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Abstract 


In this study it is shown that the delay in the spin-up of precipitation early in numerical 
atmospheric forecasts is a deficiency that can be corrected by diabatic initialization combined with 
diabatic forcing. For either to be effective requires some knowledge of the magnitude and vertical 
placement of the latent heating fields. Until recently the best source of cloud and rain water data has 
been the remotely sensed vertically integrated precipitation rate or liquid water content Still, the 
vertical placement of the condensation (heating) remains unknown. However, some information 
about the vertical distribution of the heating rates and precipitating liquid water and ice can be obtained 
from retrieval techniques that use a physical model of precipitating clouds to refine and improve the 
interpretation of the remotely sensed data. A description of this procedure and an examination of its 
three dimensional liquid water products, along with improved modeling methods that enhance or 
speed-up storm development and precipitation early in the forecast, are the subject of this report. 

The first chapter in our final report contains a description of the spin-up problem. Included is 
an examination of diabatic initialization and diabatic forcing methodologies that help enhance the spin- 
up of precipitation in numerical atmospheric forecasts. The potential usefulness of remotely retrieved 
liquid water data is also discussed. Chapters two and three of the final report describe a retrieval 
procedure that uses a physical model to help improve and interpret remotely sensed data. 

Highlighted are the physical model and the retrieval procedure, its accuracy and ability to produce 
three dimensional spatial liquid water fields using data from the Special Sensor Microwave/Imager 
(SSM/I). Retrievals from the SSM/I, qualified by information from a cloud model, are found to 
verify well against independent estimates of the liquid water. Chapter four describes an application 
in numerical model forecasts of these remotely retrieved three dimensional data sets. In the fifth and 
final chapter of the report, the governing equations for the kinematic moisture and heat flux are 
introduced, along with an examination of the influences of mesoscale diabatic processes on the large 
scale atmospheric dynamics. 
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Grant NAGW-1855 supported research for atmospheric forecast model diabatic initialization 
studies and for the development of physical retrieval methods for the Special Sensor 
Microwave/Imager (SSM/I). Funding was split evenly between these two efforts. Results from 
our research are presented in chapters 1 through 5 of this report A brief descnption/summary to 
each of the five chapters is presented next 

1. Diabatic Forcing and Initialization with Assimilation of Cloud and 
Rainwater in a Forecast Model: Methodology 

The focus of this part of the investigation is to find one or more general modeling techniques 
that will help reduce the time taken by numerical forecast models to initiate or spin-up precipitation 
processes and enhance storm intensity. If the conventional data base could explain the atmospheric 
mesoscale flow in detail, then much of our problem would be eliminated. But the data base is 
primarily synoptic scale, requiring that a solution must be sought either in nonconventional data, in 
methods to initialize mesoscale circulations, or in ways of retaining between forecasts the model 
generated mesoscale dynamics and precipitation fields. All three methods are investigated in this 
study. 

Herein we examine the initialization and assimilation of explicit cloud and rainwater 
quantities computed from conservation equations in a mesoscale regional model. The physical 
processes include condensation, evaporation, autoconversion, accretion and the removal of rainwater 
by fallout The question of how to initialize the explicit liquid water calculations in numerical 
models and how to retain information about precipitation processes during the 4-D assimilation cycle 
are important issues that will be addressed. The explicit cloud calculations have purposely been 
kept simple so that different initialization techniques can be easily and economically tested. 
Precipitation spin-up processes associated with three different types of weather phenomena are 
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examined. 


Our findings show that diabatic initialization, or a diabatic initialization in combination with 
a new diabatic forcing procedure, work effectively to enhance the spin-up of precipitation in a 
me so scale numerical weather prediction forecast Also, the retention of cloud and rain water during 
the analysis phase of the 4-D data assimilation procedure is shown to be valuable. Without detailed 
observations, the vertical placement of the diabatic heating remains a critical problem. 

\ Vt KiL J 

2. Physical Retrieval of Precipitation Water Contents from Special Sensor 
Microwave/Imager (SSM/I) Data- Part I: A Cloud Ensemble/Radiative 
Parameterization for Sensor Response 

The physical retrieval of geophysical parameters based upon remotely sensed data 
requires a sensor response model which relates the upwelling radiances that the sensor observes to 
the parameters to be retrieved. In the retrieval of precipitation water contents from satellite passive 
microwave observations, the sensor response model has two basic components. First, a description 
of the radiative transfer of microwaves through a precipitating atmosphere must be considered, 
because it is necessary to establish the physical relationship between precipitation water content and 
upwelling microwave brightness temperature. Also the spatial response of the satellite microwave 
sensor (or antenna pattern) must be included in the description of sensor response, since precipitation 
and the associated brightness temperature field can vary over a typical microwave sensor resolution 
footprint. 

In the current study a “population” of convective cells, as well as stratiform clouds, are 
simulated using a computationally-efficient multi-cylinder cloud model. Ensembles of clouds 
selected at random from the population, distributed over a 25 km x 25 km model domain, serve as 
the basis for radiative transfer calculations of upwelling brightness temperatures at the SSM/I 
frequencies. Sensor spatial response is treated explicitly by convolving the upwelling brightness 
temperature by the domain-integrated SSM/I antenna patterns. The sensor response model is 
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utilized in precipitation water content retrievals, which is die subject of chapter 3 in this scries. 

3. Physical Retrieval of Precipitation Water Content from Special Sensor 
Microwave/Imager Data- Part II: Retrieval Method and Applications 

A physical retrieval method for estimating precipitating water distributions and other 
geophysical parameters based upon measurements from the DMSP-F8 SSM/I is developed. Three 
unique features of the retrieval method are (1) sensor antenna patterns are explicitly included to 
accommodate varying channel resolution, (2) precipitation-brightness temperature relationships are 
quantified using the cloud ensemble/radiative parameterization described in chapter 3, and (3) spatial 
constraints are imposed for certain background parameters, such as humidity, which vary more 
slowly in the horizontal than the cloud and precipitation water contents. The general framework of 
the method will facilitate the incorporation of measurements from the SSM/T, SSM/T-2 and 
geostationary infrared measurements, as well as information from conventional sources (e.g., 
radiosondes) or numerical forecast model fields. 

4 . Forecast Model Applications of Retrieved Three Dimensional Liquid Water 
Fields 

Forecasts are made for tropical storm Emily using heating rates derived from the SSM/I 
physical retrievals described in chapters 2 and 3. Average values of the latent heating rates from the 
convective and stratiform cloud simulations, used in the physical retrieval, are obtained for individual 
1 . 1 km thick vertical layers. Then, the layer-mean latent heating rates are regressed against the slant 
path-integrated liquid and ice precipitation water contents to determine the best fit two parameter 
regression coefficients for each layer. The regression formulae and retrieved precipitation water 
contents are utilized to infer the vertical distribution of heating rates for forecast model applications. 
In the forecast model, diabatic temperature contributions are calculated and used in a diabatic 
initialization, or in a diabatic initialization combined with a diabatic forcing procedure. Our forecasts 


3 



show that the time needed to spin-up precipitation processes in tropical storm Emily is greatly 
accelerated through the application of the data. 

5. Moist Wind Relationships 

We have identified the equations describing the temporal and spatial behavior of the kinematic 
moisture and heat flux. In these nonlinear equations, the contribution by diabatic processes to the 
large-scale flux is composed of two parts. One part is associated with a Rayleigh damping term 
while the other arises from temporal and spatial changes in the pressure gradient term. Given 
below are some findings from the manuscript accepted for publication in Mon . Wea. Rev.. 

We have found that the influence of diabatic processes on large-scale moisture fluxes depends 
greatly on the degree of balance between forcing and damping terms in the governing kinematic flux 
equations. The existence of a near balance requires a reduction in the large-scale horizontal 
geostrophic wind speed. Based on an examination of the moisture flux equations, it is argued that 
reductions in the large-scale horizontal wind speed observed within major cumulus cloud systems 
help conserve large-scale moisture fluxes. The deviation of the wind from geostrophic conditions 
is easily estimated for a near balanced state. This wind modification induces secondary vertical 
circulations that contribute to convergence, creating or supporting long-lived mesoscale flows. We 
believe this process to be a major supporter of the mesoscale circulations observed in severe storms 
and squall lines. In the tropics the wind modification has an antitriptic relationship. 

These diagnostic findings suggest possible modifications to the wind field in the application of 
a cumulus parameterization, and may be important in diabatic initialization of numerical weather 
prediction models. 
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Abstract 


In this study diabatic initialization and assimilation techniques are tested in a semi-implicit 
hydrostatic forecast model containing explicit representations of cloud and rain water. Diabatic 
vertical mode initialization and a diabatic forcing process are combined to strengthen the diabatic 
signal and to help the forecast keep the modification. In some experiments explicit cloud and rain 
water fields are retained through the analysis/initialization phase of the forecast cycle, i.e., from the 
end of one forecast to the beginning of the next sequential forecast, to help reduce the spin-up time 
of major precipitating events. This is only possible when these variables are explicitly computed 
and retained for use in prognostic cloud and rain water calculations. Deficiencies in the forecast 
model’s production of precipitation can be reduced provided liquid cloud and rain water data are 
available. 

Numerical experiments for a wide range of grid scale sizes, varying from 150 to 40 km, 
illustrate the adaptability of the diabatic initialization, diabatic forcing and liquid water assimilation 
procedures in combination with explicit cloud and rain water calculations. Particulars include 
information about the simple explicit cloud scheme and the diabatic initialization and diabatic 
forcing procedure. Presented also are comparisons between forecasts computed with the diabatic 
initialization, diabatic forcing and retained cloud and rain water fields and those obtained by the 
traditional adiabatic approach. 
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1. Introduction 

In this investigation the primary goal is to identify one or more general procedures that will 
help reduce the time taken by numerical weather prediction models containing explicit cloud and 
rain water formulations to initiate or spin-up precipitation processes. These precipitation 
intensifying procedures also assist in enhancing the storm intensity. If the conventional data base 
could explain the atmospheric mesoscale flow in detail, then much of our problem would be 
eliminated. But the data base is primarily synoptic scale, requiring that a solution must be sought 
either in nonconventional data, in methods to initialize mesoscale circulations or by procedures that 
retain the model generated mesoscale dynamics and precipitation fields between the end of a 
forecast and the beginning of the next sequential forecast All three methods are investigated in 
this study. These methods would allow the existing nonconventional data base, including satellite 
estimates of precipitation rate, to become more valuable for fine tuning and correcting model 
inaccuracies and conventional data shortcomings. 

In this study we examine in a mesoscale regional model the initialization and assimilation of 
explicit cloud and rain water quantities. Parameterizations of the physical processes of 
condensation, evaporation, autoconversion, accretion and the removal of rain water by fallout are 
all included as source and sink terms in the liquid water conservation equations. The question of 
how to initialize the explicit liquid water calculations in numerical models and how to retain 
information about precipitation processes during the 4-D assimilation cycle are important issues 
that will be addressed. The explicit cloud model has purposely been kept simple so that different 
initialization techniques can be easily and economically tested. 

With increasing computer speed there has been a slow but steady enhancement in model 
horizontal grid resolution. Current global numerical weather prediction models are approaching 
the spatial resolution where time-dependent cloud variables are physically relevant The time to 
study the feasibility of initializing and assimilating cloud and rain water using explicit cloud 
conservation equations within a regional mesoscale model is now. The knowledge of explicitly 
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calculated cloud and rain water fields as a function of time and space provides information, such as 
model heating rates and vertical distributions, that are more difficult to obtain from conventional 
cloud parameterization approaches where cloud properties are not retained or used in a prognostic 
manner. 

Goud and rain water are the cumulative products of the transport of atmospheric moisture, 
both vapor and liquid, plus the physical processes included in cloud development and decay. By 
their very nature cumulative processes are difficult to initiate in numerical models, and it is 
common for precipitation to be absent during the first few hours of a numerical forecast One 
reason for this deficiency is the resolution of the data. Conventional surface and upper air reports 
are frequently unable to capture the essence of mesoscale circulations with their enhanced 
convergence which supports the mesoscale precipitation process. The analysis of the irrotational 
wind component in the tropics is especially poor (Mizzi and Kasahara 1989). In mid-latitudes the 
large scale part of the irrotational wind can be estimated fairly accurately (Kasahara et al. 1992) but 
the mesoscale component depends a great deal upon the diabatic term (Tarbell et al 1981). The 
missing mesoscale circulations also influence the moisture and thermal fields (e.g., Krishnamurti et 
al 1984; Turpeinen 1990). Consequently, in order to initiate precipitation and assimilate cloud 
parameters it is necessary to know additional information beyond the conventional synoptic data 
base. This might include accurate estimates of the vertical heating rates and the low level 
convergence and/or the high level divergence that produces the cloud-sustaining vertical motion 
field. These also include the supporting mesoscale circulations that help produce long-lived 
disturbances. The incorporation of this data into model simulations is further complicated by 
inaccuracies or shortcomings in the forecast, e.g., phase errors. 

When the heating rates are known or can be approximated, then either diabatic (Errico and 
Rasch 1988) or cumulus initialization (Donner 1988) is possible. These initialization procedures 
help establish circulation patterns consistent with latent heating processes that are otherwise not 
resolved in synoptic scale data. For a review of some studies using diabatic initialization schemes, 
see Table 1 in Turpeinen et al. (1990) and the introduction of Mathur et al ( 1992). Without the 
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heating rate information each forecast must go through a short time period where the mesoscale 
circulations that support precipitation processes must 'spin-up' before rainfall can be produced. In 
regional mesoscale models this spin-up usually occurs during the first 2 to 6 hours. Diabatic 
initialization should reduce the degree and duration of the spin-up, but this procedure is not always 
completely successful (Errico and Rasch 1988; Kasahara et al. 1992). In contrast, normal mode 
or dynamic initialization procedures, designed to help with the removal of gravity waves during 
initialization, only provide slight or no assistance in speeding up precipitation processes (Lejenas 
1980). 

Goud and rain water data are generally not available. Rain gauge reports imply information 
about the total amount of liquid water reaching the surface, but no knowledge about the 
instantaneous vertical distribution is recoverable. Most microwave satellite measurements of 
instantaneous cloud quantities only resolve the vertically integrated cloud and/or integrated rain 
water. Thus it is necessary to make some assumption about the vertical distribution of the heating 
rate. New retrieval methods that combine observed microwave satellite radiances with cloud 
models should provide some insight into the vertical partitioning of the liquid water. These 
techniques are just beginning to be examined (Olson 1993). 

The use of diabatic initialization requires an answer to the question of how the latent heating 
is going to be distributed in the vertical. Some studies have used representative profiles (Fiorino 
and Warner 1981; Puri and Miller 1990) others assume a general parabolic profile (Tarbell et al 
1981; Salmon and Warner 1986, etc.) while some have tried to use the forecast model itself in 
some fashion with observations (Wergen 1988; Dan aid 1985; Turpeinen et al. 1990; Donner and 
Rasch 1989; Puri and Miller 1990; Heckley 1990). Donner (1988) used a variational formulation 
to optimize the vertical distribution with respect to the methodology of the Kuo cumulus 
parameterization. In all but the Donner study, it has been customary to use either the diabatic 
omega equation or a diabatic form of the normal mode initialization scheme. In contrast, Donner 
used a variational optimization. 
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Early pioneering work by Krishnamurti et al. (1984) suggested that in addition to the diabadc 
heating a re-analysis of the moisture and wind field is needed. Even earlier studies by Perkey 
(1976) and Wolcott and Warner (1981) enhanced the moisture field in cloudy regions. To 
produce the observed rainfall, Donner (1988) made slight modifications in the temperature field 
and relatively large changes to the humidity field. Donner and Rasch (1989) found that die 
divergence, temperature and moisture must all be adjusted to produce the proper precipitation. The 
necessity of moisture initialization to reduce the spin-up problem has also been noted in Turpeinen 
(1990) and Puri and Miller (1990). 

The cumulus initialization used by Donner (1988) is effective in dealing with cumulus 
convective situations. However, improvements are still needed in the initialization of grid scale 
precipitation processes. In this study we use diabatic initialization, diabatic forcing, allow for 
moisture modification and in some experiments assimilate cloud water and rain water predicted by 
a forecast beginning 12 hours earlier. For the latter cases, non-zero cloud and rain water are 
carried through the 4-D data assimilation process from the end of one forecast to the beginning of 
the next forecast The diabatic initialization is installed to work with the vertical mode initialization 
scheme of Bourke and McGregor (1983). Diabatic forcing is formulated in a new design to 
represent a natural continuation of the diabatic initialization. Together, diabatic initialization and 
diabatic forcing form a powerful system to reduce the time needed to spin-up precipitation 
processes and enhance storm development. 

Theoretically, if precipitation data were always available then cloud and rain water could be 
treated much like other meteorological variables, i.e., updated and modified during the 4D 
assimilation process. However, caution is required when making large changes to existing model 
cloud fields, or to the dynamics that support the precipitation process. Testing and experience with 
explicit cloud initialization will dictate what the correct procedure should be. The problem goes 
far beyond any single study, and must be addressed in an operational setting where a large number 
of cases and situations can be examined over a long time period. Because of resource limitation 
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this study tests specific parts of the proposed initialization procedures separately, and then only on 
a few cases. The sensitivity of the initialization to the vertical distribution of the diabatic heating 
is appraised. Emphasized are procedures that have the potential to reduce the spin-up time and 
assist with enhancing storm development. 

2. Modeling Considerations 

A modified version of the Bureau of Meteorology Research Centre (BMRC) operational 
incremental O/I analysis (Mills and Seaman 1990), vertical mode initialization (Bourke and 
McGregor 1983), and semi-implicit forecast model (Leslie, et al. 1985) are used in this 
investigation. The incremental analysis is set to preserve the irrotational wind field while the 
rotational component of the differences between data and model guess wind and other fields are 
evaluated by the O/I. In this study the horizontal grid configuration varies from 40 to 150 km, 
while the vertical sigma coordinate is subdivided into 15 levels. The horizontal wind components, 
temperature and mixing ratio fields are smoothed with the 6th order implicit low-pass tangential 
filter (Raymond 1988), with a filter parameter setting sufficiently large enough to provide some 
horizontal mixing. Cloud and rain water quantities are not filtered, nor is horizontal diffusion 
applied to these fields. 

Explicit conservation equations for cloud water and rain water are inserted in the model. 

The explicit cloud physics contain processes modeled using a modified Kessler type scheme 
(Kessler 1974), while the condensation and fall velocity calculations are similar to those described 
in Anthes, et al. (1987) and Liu and Orville (1969). The 'fallout' term is calculated by a semi- 
Lagrangian method to allow for the larger time steps used in the semi-implicit model calculations. 
Ice phase cloud physics have not yet been explicitly included in the model, but modifications are 
incorporated to increase the life time of high clouds and and to enhance anvil structures. 

In sigma coordinates the frictionless explicit conservation equations for water vapor, cloud 
and rain water satisfy 
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(la) 

9(p*q v )/9t + 9(up*q v )/9x + 9(vp*q v )/9y + 9(op*q v )/9o= P*( p evap ” P COn) + ’ 


(lb) 

9(p*q c )/9t + 9(up*q c )/9x + 9(vp*q c )/9y + 9(6p*q c )/3c^= - p*(P AC c + p auto~ p con) 

+ p z » 


(lc) 

9(p*q r )/9t + 9(up*q r )/9x + 9(vp*q r )/9y + 9(op*q r )/9o= P*( p A cc + P AUTO P CON) 

-g9(pq r v t )/9a + F z . 


Here q v , q c and q r are the water vapor, cloud water and rain water mixing ratios. p evap * s 
evaporation of raindrops, p con the condensation of water vapor (Asai 1965) or evaporation of 
cloud droplets, p ACC is the accretion of cloud droplets by raindrops and p AUTO represents the 
autoconversion of cloud droplets into raindrops. The acceleration of gravity is g, p* is the 
surface pressure, p is the density of dry air and v t is the mean terminal velocity of the raindrops. 
The independent variables t, x, y, and c have their standard definitions. F z represents the 
vertical diffusion. The map factor is omitted in the above equations. 

Because of horizontal scale considerations, our explicit cloud calculations do not model 
individual clouds but yield approximate average values for a system or complex of clouds. This 
was the motivation behind the modifications made to the cloud physics formulation. It is known 
that the lifetime of individual convective clouds is less than one hour while a cloud complex exists 
for many hours with remnants lasting for days. Thus Kessler's (1974) accretion (collection) 
process was modified to slow the conversion from cloud to rain water when the cloud water q c is 
below a preset critical value qc^^r Also incorporated into the accretion process is a vertical 
dependence on temperature that follows the scheme suggested by Schlesinger (1990). These 
changes help retain and enhance anvil or cloud top features. In addition to explicitly modeling 
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cloud complexes, the consequences of convection are parameterized by a modified Kuo-type 
cumulus parameterization scheme (Sundqvist et al 1989). However, in some cases the cumulus 
parameterization scheme is intentionally turned off to assist in the study of the spin-up of grid scale 
precipitation. 

In mixing ratio units (gg 1 ), the formula for the rate of change from accretion is 


Pacc^^C^ 9c 9 r 


0.875 


( 2 ) 


where 


C^minfl.O, 1.0 - (C^, - %)/ C CTit ] 

C2=max[0.0, min(1.0, (T - 233)/ 40)] . 

Here T is the temperature in degrees Kelvin, and C cri =l gkg 1 . Comparisons between model 
cloud water fields and infrared VISSR Atmospheric Sounder (VAS) images of cloud patterns and 
categories of low or high clouds, and SSM/I retrieved precipitation rates indicate acceptable ranges 
of cloud and rain water predicted by the explicit formulation. This will be shown later. 

The autoconversion process is described by 

^ > auto = ^i( 9c ~ 9co) » 9c - 9co , 

(3) 

Pauto^ > 9c < 9co • 

The parameters k^lO' 3 s' 1 and the onset of autoconversion depends on the critical value which 
has the value of 0.2 g kg' 1 . 
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To facilitate the large time step used in the semi-implicit forecast model die fallout term in (lc) 
is rewritten as two terms, one of which represents vertical advection, i. e., 

g9(pv t q r )/a<j a Aa(p*q r )/3o + (p*q r g/p*)3(pv t )/3o . (4) 

Here A (s' 1 ) is equivalent to a velocity in the vertical a coordinate system, that is 

A= gpv/p* . (5) 

Thus, the time tendency for p*q r in (lc) is found by semi-Lagrangian procedures using (4) and (5). 

3. The Diabatic Initialization and Diabatic Forcing Formulation 

The diabatic initialization is installed into the forecast model to take full advantage of the 
semi-implicit method. This is accomplished by modifying the semi-implicit representation of the 
average temperature T 2t , computed from the thermodynamic equation (McGregor et al 1978; Leslie 
et al 1985), by adding the temperature contribution T^ from a measured or approximated diabatic 
heat source. This redefines the average temperature at the appropriate grid points using 

T 2t = T 2t +Tex (6) 

(This is approximately equivalent to adding to the thermodynamic equation a forcing term T EX /At, 
and including it in the semi-implicit formulation.) This modified temperature is then used in the 
semi-implicit calculations of the finite difference averaged horizontal wind components and surface 
pressure. To avoid excessive evaporation, some upward adjustment of the relative humidity may 
be required during this process . The amount of moisture needed to satisfy the expected 
condensation requirements is c p T EX /L, where c p is the specific heat for constant pressure and L is 


10 



the latent heat of condensation. When cloud and rain water are assimilated it is especially 
important to retain saturation at the appropriate levels. During the vertical mode initialization 
procedure (Bourke and McGregor 1983), Tgx is held constant at the selected grid points. Recall, 
that the averaging in (6) for the semi-implicit method, with any variable x» is defined by 

X 2t = [x t + i + x x-i]/2 . (7) 

Following the initialization, the immediate application of an external diabatic forcing 
represents a natural extension of the diabatic initialization. Wang and Warner (1988) showed that 
this was useful. To keep the forcing from becoming excessive, the grid point average temperature 
is modified by the difference between the external contribution Tgx and the model produced or 
internal diabatic temperature change T^AtLQ/Cp, yielding 

T 2t sT 2t + (Tex-Tn). (8) 

Here Q is the model produced condensed rate and At the time step. In our study the quantity in 
parenthesis in (8) is taken to be non-negative. If exceeds % at a given grid location and 
timestep, then the forcing is reduced to zero. At some time in the forecast this would be the 
desired result but, it is not explicitly required before the termination of the forcing at a preselected 
forecast time. Too much forcing can over-develop the updraft and the liquid water content and/or 
create excessive gravity waves. The inclusion of the internally generated latent heating in (8) 
helps reduce over-development. 

Note that (6) and (8) are similar, except that in (8) model produced condensation reduces the 
external forcing. This diabatic forcing is a natural continuation of the diabatic contributions used 
in the initialization. Thus, it is applied starting with the first time step after the vertical mode 
initialization and continued until the model produces an equivalent temperature change or until the 
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termination time is exceeded. This procedure is illustrated in Fig 1 a. 

To simulate the movement associated with a disturbance, Tgx is advected by the model 

horizontal winds according to 

3(Tex)/^ 1 + u3(T ex )/^ x + v3(TEx)/3y = 0- ^ 

In regions with strong vertical shear it may be advantageous to advect Tgx using a vertically 
averaged horizontal wind. The advection process will ultimately smooth and diffuse the heating if 
(9) is used over a sufficiently large time period. Normally, the diabatic forcing is not applied long 
enough for this to become a problem. Albeit, to avoid noise some tapering is required to terminate 
the diabatic forcing, provided (8) is not reduced to zero by model generated condensation 
processes. Remotely sensed time variations in Tgx could also be specified if these were known. 

In their absence, the advection of T^ is necessary to approximately locate the condensation. In 
our calculations T EX remains unchanged except for horizontal advection. Thus, in this study the 
application of the external forcing is somewhat different from the procedure used in Newtonian 
nudging (Hoke and Anthes 1976). 

An alternative approach, examined in this study, is to apply the external diabatic forcing at the 
satellite observation time without reinitializing the model, i.e., with no diabatic initialization. Then, 
during the forecast, a gradual transition into and exit from the external forcing greatly avoids the 
production of noise and gravity waves. This is shown graphically in Fig lb to contrast how this 
procedure differs from diabatic initialization. The time span needed to produce gradual changes 
depends upon the magnitude of the external forcing. The length of the time period used in the 
forcing can be varied to achieve the desired result. Excessively long periods of forcing should be 
avoided since they can can lead to some negative consequences. For example, overproduction of 
precipitating liquid in the atmospheric column leads to down bursts accompanied by low-level 
divergence, ultimately choking off the moisture source. 

In any event, the horizontal distribution of the forcing is determined using the gnd 
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configuration and knowledge of the location of the storm or liquid precipitating water. In die 
absence of measurements, the best the forcing remains unknown and is a 

function of several factors, including the nature of the disturbance. For example, hurricanes are 
known to be warm core; warming extends to surface, while most middle-latitude disturbances have 
upper level lows with warming that does not extend to the surface (Wallace and Hobbs 1977). 

Raymond (1993) has identified diabatic wind modifications that induce secondary vertical 
circulations that contribute to the convergence, helping to create or support long-lived mesoscale 
flows. An approximate representation for the secondary flow within the diabatically active part of 
the disturbance, based on near conservation of the large scale kinematic moisture flux, requires that 
the x and y coordinate horizontal wind components u and v be modified by the addition of a 
perturbation term of the form 


u pert =R l V g- R l 2u g + 


W = ~ R l U g ~ R l 2v g + 


( 10 ) 


(ID 


Here Rj=— lQ/(fq) I, Q is the condensation rate, f is the Coriolis parameter, q the mixing ratio and 
u g and v g are the geostrophic wind components for u and v respectively. Generally Rj is small, 
say, 0<R 1 <-0.50. Because (10) and (1 1) only describe the modifications within the updraft, 
they represent only part of an unbalanced ageostrophic approximation of the mesoscale circulation. 
Theoretically, a balanced and complete flow can be determined, but this requires additional 
computational resources to solve a vorticity equation. In our calculations, Rj in (10) and (1 1) is 
used as an independent gauge to measure the production of latent heating, but the wind field is not 
modified. 
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4. Precipitation Events 

Precipitation spin-up processes associated with three different weather phenomena are 
examined. The first situation includes a tropical disturbance known as Emily. The initial 
conditions are taken from NMC’s initialized global data set beginning at OOOO UTC 21 Sept. 
1987. At that time Emily was a tropical storm located among the Windward Islands a little ways 
south of the Island of Martinique, just north of the Venezuelan coast in the Atlantic Ocean. The 
model calculations are performed on a 65 x 65 horizontal grid with 40 km grid spacing, having 15 
vertical sigma levels. The time step size used in the semi-implicit model is 300 s. 

The second major precipitation study examines rapid cyclonic development in the North 
Atlantic. Intensive observation period number five (IOP-5) is selected from the Experiment on 
Rapidly Intensifying Cyclones over the Atlantic (ERICA). It represents a case where a major 
cyclone over water experiences a rapid pressure deepening. Our examination extends from 0000 
UTC 19 January 1989 through 0000 UTC on the 20th. For this part of the study the forecast 
model uses a horizontal grid resolution of 120 km. The initial meteorological fields for our 
ERICA IOP-5 case study are obtained from NMC's initialized global analyses. Conventional data 
are limited because the disturbance is located over water. However, observations during ERICA 
have documented the cyclone's path, pressure and provided an assortment of other atmospheric 
measurements using dropsondes and aircraft. Also, the microwave measurements of the 
precipitation rates retrieved from the Special Sensor Microwave/Imager (SSM/I) provide a product 
of major importance in verifying the precipitating liquid in the forecast. 

The final precipitation event used to study model spin-up comes from a time period that 
includes frontal development associated with strong wind jets over the US. Specifically, 24h 
forecasts are made for the period beginning 1200 UTC 21 April 1988. This case over North 
America was selected because of its large conventional data bases including surface, upper air and 
aircraft winds. The horizontal grid scale used in this part of the modelling effort is 150 km on a 
35 x 35 grid with 15 vertical sigma levels. 
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5. Discussion of Results 

A. Tropical Storm Emily 

Super-imposed over the validating GOES infrared (IR) image of tropical storm Emily in Figs 
2 and 3 are the 700 mb geopotential height and divergence fields, respectively, from a fourteen 
hour forecast beginning at 1000 UTC 21 September 1987. The forecast uses an adiabatic 
initialization with four vertical modes and the initial meteorological fields come from a ten hour 
preforecast. Decreasing the number of modes by one or two makes no significant difference in the 
adiabatic initialization. Note in Figs 2 and 3 that the height and divergence fields show very little 
development whatsoever in the immediate vicinity of Emily. A small trough is located to the east 
of Emily in Fig 2, but no precipitation is produced over any water surface, although moderate 
amounts of precipitation fall in some elevated regions. 

In the absence of the model’s inability to produce any diabatic heating in the vicinity of Emily 
during a 10 h preforecast beginning 0000 UTC 21 Sept 1987, a proxy latent heating field is 
assigned based on the approximate position of Emily in the 1001 UTC satellite image. (In a 
separate study that will be presented elsewhere, data from SSM/I estimates of the precipitating 
liquid and ice at 1000 UTC 21 Sept 1987 are used in initialization experiments) In this study 
highly simplistic initial conditions are used to test our basic diabatic forcing procedure. Previous 
investigators have also used simplified approaches to describe the vertical distribution of the latent 
heating, including Tarbell et al (1981), Salmon and Warner (1986), Ninomiya and Kurihara (1987), 
and Wang and Warner (1988). Our proxy temperature distribution is described by the product of 
vertical and horizontal sine functions, viz., 

t ex =K sin(a7t)sin(b7t)sin(c7t) (12) 

where K is the maximum magnitude of T E x and a=(k-k T )/(k B -k T ), b=(j-j s )/(jN-j s ) and 
c=(i-i E )/(i w -jE)- Th e independent coordinates i, j, and k are associated with the x, y and a 
variables, respectively. At 1001 UTC 21 September 1987, Emily was approximately centered 
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near grid point location (40.5,30.5) in our 65 x 65 horizontal grid. Assigning jj— 26, j^— 35, 
i E =37, i w =44 and setting K=0.5 gives the distribution shown at 500 mb in Fig 4, provided the 
sigma levels are between 0.25 and 0.80. This vertical distribution is compatible with pre- 
hurricane/tropical wave observations (Kuo and Raymond 1980; Shapiro and Stevens 1980). In 
general, the values for k B and k T are varied to test the impact of heating at different vertical heights. 

Figures 5 through 7 contain fields from a fourteen hour forecast using the diabatic 
initialization and a two hour application of diabatic forcing as described in (6) through (12), and 
illustrated in Fig la. For convenience, ^iabaticjnitialization will be denoted by DI, while DF(K) 
will identify ^iabaticjforcing with a maximum temperature perturbation of K degrees. The 700 mb 
height field from the fourteen hour forecast, Fig 5, contains closed contours of small magnitude 
over the southern extent of tropical storm Emily, indicating that the forecasted disturbance still has a 
cold core. The 700 mb vorticity, with maximum in excess of2x 10 4 s 1 is similarly located (Fig 
6), while the 700 mb divergence (Fig 7) is positive directly behind the disturbance and negative 
(convergence) within the vicinity of the storm and to the southwest The diabatic temperature 
perturbation used in these calculations has a maximum value of K=0.5 C per time step of 300 s. 

Six iterations (iter=6) are used in the vertical mode initialization (Bourke and McGregor 1983), each 
having a time step of 50s, requiring that K be modified appropriately in the initialization, i.e., 

K=0. 5/iter. In (12), the parameters describing the vertical depth of the perturbation are assigned 
=9 and kj=3, i.e., the forcing occurs between sigma levels 0.25 and 0.80. Also, the diabatic 
forcing, which is applied over a time period of two hours, is ramped linearly to zero during the 
second hour. 

Four vertical wave numbers are used in the initialization shown in Fig. 8. Included in this 
figure are curves describing the number of grid scale condensation points per time step for 
adiabatic, diabatic initialization (DI), and for two different time periods of diabatic forcing (DF). 
Note that, without the diabatic forcing, the diabatic initialization is unable to produce sustained 
condensation activity over the fourteen hour forecast This is true for vertical modes 2 and 3 also. 
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Adding additional moisture during the initialization at all locations influenced by the diabatic 
heating does not enhance the long term development. However, strong diabatic activity occurs 
when the diabatic initialization is combined with diabatic forcing, where the forcing is applied for a 
duration of either two or three hours. In Fig 8 note that the longer the diabatic forcing is applied, 
the more intense is the storm development as gauged by the production of precipitation. Both the 
diabatic initialization and diabatic forcing are restricted to the immediate Emily region which is 
entirely over water, as illustrated in Fig 4. Thus, the difference between the adiabatic and other 
curves represent development in the Emily region. In fact, the development wants to become too 
extreme, partially because of the large availability of moist air in this tropical location, and partially 
because the natural tropical storm circulation, feedback and balancing processes, developed over 
several days, are difficult for the model to duplicate in a few hours. To keep the storm from over 
developing or developing too rapidly, the maximum amount condensed at any one time step at any 
grid location is restricted by requiring that -0.5<Ri<0. Conservation of moisture is retained since 
the excess is carried in the slightly super saturated mixing ratio. This is thought to be a better 
alternative than applying artificial diffusion which is non-conservative. It is also consistent with 
the fact that small sub-grid scale moist and dry downdrafts exist within the cloud environment, 
limiting to some extent the updraft while enhancing mixing, thereby reducing the condensation 
amount. Except for the K-theory vertical diffusion, these processes are not explicitly included in 
the conservation equations (la,b,c) or in the condensation calculation, but should be parameterized 
in some way. 

In our calculations the horizontal wind is required to spin-up through time since the 
disturbance is initially absent, and since no bogus tropical storm wind field is inserted during the 
initialization. At non-tropical latitudes it has been shown that a perturbation spiral feature 
associated with a mean converging flow has an e-folding times in excess of 22 hours (Raymond 
1989) as a consequence of vortex stretching (mean divergence of 10" 5 s). Still, in model 
simulations performed by Kurihara and Tuleya (1981), the development of a deep vortex 
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characteristic of a tropical storm from a shallow easterly wave required three to five days. 

The mesoscale circulation pattern associated with the tropical storm also influences the 
diabatic activity. When the relative vorticity is of the same order as the planetary vorticity, then the 
surface pressure tendencies are influenced very little by diabatic heating (Delden 1989). But, the 
more intense the vortex, the greater the response to a given heat source (Schubert and Hack 1982; 
Delden 1989). 

Figure 9 shows that when diabatic forcing is used, the number of vertical modes in the 
initialization has little influence on the grid scale condensation. But the magnitude of the forcing is 
shown to be influential, as illustrated by the difference between the DF(0.3) and DF(0.5) curves. 
Clearly, the greater the forcing the greater the response. Also in Fig 9 it is seen that when the 
forcing is closer to the surface, i.e., between sigma levels 0.4 and 0.95, the response is greater, 
presumedly because more moisture is transported vertically. 

In Fig 10 a maximum of 14.26 cm of precipitation is recorded over fourteen hours when the 
forcing DF(0.5), centered between sigma levels 0.25 and 0.80, is applied for two hours. 
Approximately fifty percent more rainfall is recorded when the forcing is extended in duration from 
two to three hours. As expected, reducing the maximum magnitude of the forcing from K=0.5 C 
per 300 s time step to 0.3 C reduces the rainfall, yielding a maximum of 4.59 cm for DF(0.3). 
When the forcing is placed lower in the atmosphere, i.e., between sigma levels 0.4 and 0.95, then 
more moisture condenses and more total precipitation falls, but for this case the maximum rainfall 
is nevertheless reduced. Generally, Fig 10 shows that a high degree of correlation exists between 
the number of grid scale condensation points, shown in Figs 7 and 8, and the maximum rainfall. 
To help visualize the precipitation pattern the rainfall recorded on the surface, associated with a two 
hour application of DF(0.5), is shown in Fig 1 1 . The convective rainfall resulting from the 
cumulus parameterization is small in this region (not shown), although the 40 km grid resolution is 
somewhat boarder line for the effective use of a Kuo-type cumulus parameterization scheme. 

The diabatic forcing can also be used at any time during a forecast as shown graphically in 
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Fig lb. For example, results from an application of diabatic forcing beginning ten hours into a 
twenty-four hour forecast, as the sole means to enhance the spin-up of tropical storm Emily, are 
illustrated in Figs 12 through 15. Beginning the application of diabatic forcing DF(0.25) in 
nearly the middle of a forecast at the satellite observation time, with the forcing centered between 
0.25 < a< 0.80, requires a gradual approach into and decline out of the forcing regime to avoid 
introducing noise. A time period of one and one-half hours (time steps 120 to 138) is selected for 
linearly introducing the forcing, which remains constant for one and one quarter hours (time steps 
138 to 153), while one hour (time steps 153 to 165) is used to turn off the diabatic forcing. The 
resulting geopotential heights contain closed contours (Fig 12), while significant vorticity (>2.5 x 
10 4 ) and convergence are revealed in Figs 13 and 14 respectively. Convergence even exists in 
the region containing the spiral arm in the satellite image. The placement of these features is fairly 
good. A maximum rainfall of 94.5 mm is indicated in Fig 15. In the control, no rainfall is 
produced in the region containing Emily, but over 90 mm falls on the mountain slopes. Figure 15 
also shows that doubling the forcing doubles the rainfall while low-level forcing gives a large 
maximum of 277.2 mm. The areal distribution of the precipitation produced by the twenty four 
hour forecast, using DF(0.25) (solid line), and the low level forcing (dashed line), is given in Fig 
16 and are somewhat reduced in comparison with the precipitation from the fourteen hour forecast 
shown in Fig 11. 

The number of grid scale condensation points as a function of time is shown in Fig 17 for the 
control (no forcing), and for three different vertical distributions and two different magnitudes of 
diabatic forcings. The grid scale condensation points occurring in the control case are associated 
with convection over elevated topography, none occurs in the vicinity of Emily. Differences then 
between the control and forced cases are associated with Emily. Figure 17 clearly shows that the 
magnitude of the forcing is important. 

Presented in Fig. 18 is the integrated precipitating water, in kgnr 2 , predicted by the fourteen 
hour forecast It can be compared with the SSM/I retrieved estimates (Olson 1993) given in Fig 
19. Note the latter is distributed over a much larger area and contains a comma shaped pattern, 
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while the forecasted amount has concentric contours much like the shape of the forcing. The 
maximum produced by the forecast is slightly too far south. Enlarging the initial areal size of the 
diabatic heating results in a proportional increase in the integrated rain water in the fourteen hour 
forecast. When the initial conditions are shifted to the north, then the development in the forecast 
responds in a corresponding manner (not shown). 

B. ERICA IOP-5 

This case differs from the Emily situation in that the 12 hour preforecast produces diabatic 
activity over the region of interest. However, diabatic heating fields from the last time step of the 
12 hour preforecast are thought to be unrepresentative enough that they should not be used directly 
in the diabatic initialization and diabatic forcing formulation. Consequently, the T^ fields, with a 
maximum of <1 C during a 600 s time step, are set proportional to the coverage and vertical 
distribution of the cloud water fields (not shown) predicted by the 12h preforecast. It is assumed 
that the distribution of cloud water is more representative of the distribution of heating than the 
larger size rain water particles which fall out as precipitation. 

Figure 20 illustrates, for a grid resolution of 120 km at 0000 UTC 20 1988, that the diabatic 
initialization produces a doubling of the negative omega values (solid lines) found at 700 and 400 
mb in the global analysis. This initialization did enhance slightly the storm development and 
rainfall distribution in the 12h forecast, but the predicted 988 mb low pressure is far from the 
observed 968 mb measured in this rapidly deepening storm (IOP-5), as indicated in the ERICA 
Field Phase Summary (Hartnett et al. 1989). Thus, given the model’s 120 km horizontal 
resolution, additional measures are required to better forecast this intensely developing storm. 

Results from a 24h forecast, ending at 1200 UTC 20 Jan 1988, are compared in Figs 21, 22 
and 23 against a control twelve hour forecast and with two diabatically forced cases. The diabatic 
forcing is applied to the global initial conditions, and in all calculations the cumulus 
parameterization scheme has been turned off to isolate the spin-up process for the grid scale 
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calculations. In Fig 21 the wind speed at 500 mb is shown. Note that the jet develops a 
maximum speed of 40 + ms 1 in all cases except the lower right figure, which is from a forecast 
using four hours of diabatic forcing. This forecast produces 50 + ms' 1 wind speeds in the dry slot 
just behind the upper part of the comma. Note that two hours of diabatic forcing did increase the 
area under the 40 ms 1 curve (lower left). The change in the wind field can be directly associated 
with the enhancement in the comma shape of the mixing ratio shown at the 700 mb level in Fig 22. 
Note also in the two and four hour diabatically forced forecasts, that the gradients of mixing ratio 
are largest to the rear of the storm. 

Fig 23 shows the precipitation (mm) accumulated during a twelve hour period. The final 
twelve hours of the twenty four hour forecast produced a small region with a maximum in excess 
of 70 mm, while the diabatically forced forecasts show a slight reduction in the maximum 
precipitation from 30 + mm (upper right) to 20 + mm. This small and temporary reduction is 
thought to occur because the continual diabatic forcing builds a deeper system containing more 
precipitating liquid aloft. If, however, the forcing is too strong or applied too long, then the 
upward vertical motion becomes large, causing an excessive amount of water to be suspended in 
the atmospheric column. Then, the eventual fallout of this rain water can lead to negative 
consequences, and possibly the premature death of the storm. 

The area averaged precipitation rate as a function of time is presented in Fig 24 for five 
forecasts and from satellite measurements. The latter is computed using a technique (Martin et al 
1990) that estimates the hourly precipitation rate from satellite infrared image cloud-top 
temperatures by a regression power-law formula. The coefficients in the power-law formula are 
obtained using maximum-likelihood correlations of relationships between cloud-top temperatures 
and a large number of raingauge precipitation rate observations (B. Goodman, personal 
communication). The variations seen in Fig 24 for the observations are small, but the extra time 
needed to spin-up the precipitation in the model is clear. In the legend in Fig. 24 “large scale” 
refers to the traditional method of instantaneously removing from the atmosphere water vapor in 
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excess of saturation. Clearly, this procedure needs less time to produce precipitation than the grid 
scale explicit cloud category which requires time for the different cloud processes to react and for 
the rain water to reach the earth’s surface. Retaining cloud and rain water unaltered through an 
assimilation cycle (heavy solid line) allows the precipitation to continue at the same rate as it had 
previously, denoted by “hot start” in the legend. Note that all methods come close to converging to 
the average precipitation rate seen in the satellite measurements. The forecast that uses diabatic 
forcing for four hours (not shown) produces just slightly less rainfall than the explicit category 
with no diabatic forcing (open circle). 

Super-imposed over the GOES 7 IR image at 1201 UTC 20 Jan. 1989 in Fig. 25 is the 
surface pressure field. This field is from a 12 h forecast that uses diabatic initialization (DI), 
consequently it deepens to only 988 mb instead of the observed 968 mb. In contrast, using the 
diabatic initialization (DI) and diabatic forcing (DF) over a 4 hour period, where during the last two 
hours the forcing is linearly tapered to zero, produces a 968 mb contour as shown in Fig 26. 
However, some gravity wave noise is also produced as can be seen by the waviness of the 
contours near the US coast and elsewhere. 

C. US Forecasts (1200 UTC 21-22 April 1988) 

The final case involves a meteorological situation where surface and upper air reports are 
plentiful. This situation was chosen to test how conventional data effects the assimilation of cloud 
and rain water quantities. Rainfall and other meteorological information for this time period is 
available in the Daily Weather Maps 1 . Generally, a stationary front is located over Colorado, 
Wyoming, and Montana, with the center of the surface low over Kansas. Another stationary front 
lying directly east of the surface low. The heaviest concentration of clouds are found to the north 
of the front. Rainfall amounts are generally light to moderate, with less than one-half of an inch 
over the last twenty-four hours. Upper level moisture and a band of high clouds have moved 

1 Daily Weather Maps, Weekly Series April 18-24, 1988, National Oceanic and Atmospheric 
Administration, National Weather Service, National Meteorological Center Climate Analysis 
Center. 
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across Mexico and into the southern US. Only a small amount of precipitation is associated with 
this upper level system. The integrated cloud water (x 100 kg nr 2 ) from our 12 hour forecast is 
shown in Fig 27, valid for 1200 UTC 22 April 1988. The predicted cloud cover is surprising 
good for the region north of the front, given the 150 km horizontal resolution use in the forecast. 
Boundary conditions influence the cloud development along the eastern, western and southern 
boundaries. The high level clouds along the southern part of the US are not predicted, partially 
because the southern boundary condition in our regional model calculations did not supply enough 
moisture. Initial conditions for this forecast used surface and upper air observations, along with 
the traditional model produced guess fields plus cloud and rain water information supplied by the 
12 hour preforecast. These liquid water quantities helped reduce the time needed for precipitation 
to develop, as illustrated below. 

Figure 28 illustrates the sensitivity of the vertical mode initialization procedure. Note that for 
two vertical modes the number of condensation points, in the grid scale calculations for a one time 
step forecast, change little when the magnitude of the heating is increased by a factor of four. 
However, the number of condensation points increases dramatically as the number of modes in the 
vertical mode initialization is increased. This requires that care must be exercised in selecting the 
proper number of vertical modes. When diabatic forcing is combined with the diabatic 
initialization, the importance of selecting the proper vertical modes is less important, as has been 
illustrated earlier. In the above calculations the Tgx fields are generated using cloud water 
information from the 12 hour preforecast. 

Changes in the grid averaged precipitation rate resulting from the inclusion of data inserted 
during the assimilation cycle are contrasted in Fig 29 with the contribution from the cumulus 
parameterization process. Here MG stands for the preforecast model generated 150 km resolution 
fields used in the initial conditions, and represents a forecast that does not use cumulus 
parameterization unless so stated. Note in Fig 29 that the largest difference in precipitation 
production between the MG simulation with cumulus (thin line) and without cumulus (circle) exists 
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early in the forecast. Also clear in Fig 29 is the fact that precipitation contributions added by the 
cumulus parameterization are larger than the precipitation changes induced by the addition of the 
surface, upper air and air craft data (thick line). To really understand the impact of data, i.e., 
temperatures, pressures, mixing ratios, and winds upon the assimilation of cloud and rain water 
will require an examination of a large number of cases. Such a study must be conducted within an 
operational environment 

The number of grid scale condensation points (Fig 30) and the grid averaged precipitation 
rate (Fig 31) indicate large differences which exist between the global no data + DI and the other 
MG categories. For this particular case, at the 1 50 km resolution, the inclusion of cloud, cloud 
and rain, or cloud and rain water plus diabatic initialization (DI) makes a difference during the first 
half of the 12 hour forecast. These differences are still much less than those found between the 
global and the higher resolution MG data sets. Nevertheless, even in this situation which 
contains weak to moderate precipitation events, there are some positive benefits, seen especially in 
Fig 31, from including the cloud water in the assimilation cycle. Their inclusion also temporarily 
increases the total moisture supply, until it all is rained out Also, diabatic initialization is usually 
required in the assimilation of cloud and rain water to help sustain or add to the supporting 
mesoscale circulation. 

The amount of mesoscale resolution in the initial model fields makes a difference 12 hours 
later in the vertical distribution of the area averaged cloud and rain water, as shown in Fig 32 and 
33. Note that the MG higher resolution initial fields produce cloud water in the 12 hour forecast 
that has a larger percentage of its total contribution at elevated sigma levels (Fig 32). Similar 
comments apply to the rain water displayed in Fig 33. This would seem to be a natural 
consequence of having enhanced vertical motions and larger magnitudes of convergence. The 
differences illustrated in Fig 32 and 33 also imply that diabatic temperature initialization, used with 
the global data is not able to correct for missing diabatic mesoscale circulation deficiencies. 
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VI Summary 

In this study several methods are shown to enhance the development and reduce the spin-up 
of precipitation when explicit cloud and rain water calculations are used. These include diabatic 
initialization, diabatic forcing and the assimilation of cloud and rain water fields from one forecast 
to the next sequential forecast To statistically validate any of the precipitation enhancing 
procedures requires an operational setting where a large number of cases can be tested. 

Diabatic initialization can be given an assist by a short period of diabatic forcing following the 
initialization. This is also beneficial for enhancing storm intensity. The forcing procedure can be 
constructed to be a natural continuation of the diabatic initialization, but can be applied at any time 
provided care is taken to avoid introducing noise. Advection of the external heating field helps 
with the proper horizontal placement of the forcing provided the forecast wind field is 
approximately correct. Extra moisture can also be added to help encourage diabatic development 
When explicit cloud schemes are employed, the assimilation of cloud and rain water from one 
forecast to the next sequential forecast has the potential to be a useful method for encouraging the 
early development of precipitation, especially in major disturbances. Improved modeling 
techniques using diabatic forcing and cloud and rain water assimilation are thus able to help reduce 
the spin-up of grid-scale precipitation and enhance mesoscale circulations, which in turn provide 
the support necessary to sustain continued diabatic activity. 

The vertical placement of the diabatic forcing remains a difficult problem. Using the vertical 
levels and proportions of the cloud development in the preforecast provides some vertical 
information, but difficulties remain in regions that lack cloud development. Most of the time the 
lack of cloud development can be traced to the lack of information in initial conditions. This study 
shows that data assimilation, diabatic initialization and diabatic forcing procedures can be used with 
explicit cloud and rain water formulations to encourage, enhance and modify the cloud and rain 
water fields, and to intensify the forecasted storm development. A major problem is the lack of 
data to initialize the cloud and rain water fields and an incomplete knowledge of the correct vertical 
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profile of latent heating. Promising satellite retrieval techniques, that combine standard remote 
sensing procedures with physical models to estimate the precipitating liquid water and ice 
distributions, may also be able to estimate the vertical placement of the latent heating fields. This 
holds great promise for providing data for the assimilation of the initial cloud and rain water fields, 
especially over water surfaces. 
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Fig. 1. A schematic showing in (a) that the diabatic forcing follows immediately after the diabatic 
initialization, while in (b) the diabatic forcing is applied later in the forecast. 

Fig. 2. Presented are the 700 mb heights from the fourteen hour forecast in the absence of diabatic 
forcing or diabatic initialization. 

Fig. 3. Same as Fig 2 but for divergence field. 

Fig. 4 Initial temperature perturbation field at 500 mb for % as obtained from (12) with K=0.5C. 

Fig. 5 Height field at 700 mb from a fourteen hour forecast containing both diabatic initialization 
and diabatic forcing DF(0.5). Contour interval is every 4 m. 

Fig. 6 Same as Fig. 5 but for vorticity. Contour interval is 4.0 x 10 5 s 1 . 

Fig. 7 Same as Fig. 5 but for divergence. Contour interval is 1.0 x 10 -5 s 1 . 

Fig. 8 The number of grid scale condensations as a function of time step is illustrated for adiabatic 
and diabatic initialization with and without diabatic forcing. 

Fig. 9 Same as Fig. 8 but showing the variation as a function of the number of vertical modes, 
vertical position and magnitude. 

Fig. 10 Bar graph of maximum rainfall for a select number of simulations. 

Fig. 1 1 Total rainfall (cm) from the fourteen hour forecast with diabatic forcing. The maximum 
amount is in excess of 14 cm. Contour interval is 2 cm. 
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Fig. 12 Height field from the twenty four hour forecast with diabatic forcing DF(0.25) inserted 

gently at time step 120 (hour 10) and removed gently by step 165 (hour 13.75). Contour 
interval is 4 m. 

Fig. 13. Same as Fig. 12 but for vorticity. Contour interval is 4.0 x 10' 5 s 1 . 

Fig. 14 Same as Fig. 12 but for divergence. Contour interval is 1.0 x 10- 5 s 1 . 

Fig. 15 Bar graph of maximum rainfall as a function of the diabatic forcing. Control case contains 
no rainfall in the vicinity of Emily but does have significant rainfall over elevated regions. 

Fig. 16 Total rainfall is shown for two cases having different diabatic forcing in the vertical. 

Dashed lines are associated with lower level forcing while the elevated forcing gives the 
solid line. Contour interval is 4 cm. 

Fig. 17 The number of grid scale condensations are shown as a function of the diabatic forcing 
which varies in its vertical distribution or magnitude. 

Fig. 1 8 Shown is the integrated rain water (kg nr 2 x 10) from the fourteen hour forecast. 

Fig. 19 Retrieved estimates of the integrated precipitating liquid and ice from SSM/I radiances 
(Olson, 1993). 

Fig. 20 Presented are the omega distributions in the adiabatic (top) and diabatic (bottom) initialized 
fields associated with ERICA IOP-5. Left-hand panels show the 700 mb distributions, 
while the right-hand panels give the 400 mb distribution. 

Fig. 21 Presented are the 500 mb wind speed from the twenty four hour (control) forecast (upper 
left), twelve hour forecast with adiabatic initial conditions (upper right), twelve hour 
forecast with 2 hours of diabatic forcing (lower left), and twelve hour forecast with 4 hours 
of diabatic forcing (lower right). 
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Fig. 22 Same as in Fig 21 but for the 700 mb mixing ratio. 

Fig. 23 Same as in Fig 21 but showing the rainfall. 

Fig. 24 Shown is the grid averaged precipitation rate (mm h 1 ) during a twenty four hour period 
ending at 1200 UTC 20 Jan 1988. The grid scale calculations with no clouds are labeled 
large scale while explicit cloud calculations, with and without a diabatic initialization and 
initialized cloud and rain water, can be compared against estimates of the hourly 
precipitation rate obtained from satellite infrared image cloud-top temperatures. 

Fig. 25 Mean sea level pressure from a 12 h forecast beginning with a diabatic initialization. 
Contour interval is 4 mb. 

Fig. 26 Mean sea level from a 12 h forecast that uses diabatic initialization and diabatic forcing. 
Contour interval is 4 mb. 

Fig. 27 The vertically integrated cloud water (x 100 kg nr 2 ) at 1200 UTC 22 April 1988 from a 12 
hour forecast initialized with cloud and rain water from the 12 hour preforecast. Contour 
interval is 3 . 

Fig. 28 The response in the grid scale condensations, as a function of vertical modes, to changes in 
the magnitude of the diabatic initialization. 

Fig. 29 Shown is the average precipitation per time step, using in the assimilation the model 

generated (MG) fields from the preforecast, with and without the cumulus parameterization 
scheme and the surface and upper air data. 

Fig. 30 The number of grid scale condensations are shown as a function of cloud water (Cl), rain 
water (Rn), diabatic initialization (DI) and whether the model generated fields (MG) or 
global fields are used in the assimilation process. 

Fig. 3 1 Same as Fig. 30 but showing the average precipitation rate. 
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Fig. 32 The per cent of the total cloud water at each vertical level at the end of a 12 h forecast is 
shown for two forecasts that have different amounts of data in the initialization. 


Fig. 33 Same as Fig. 32 but for total rain water. 
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1. INTRODUCTION 

The physical retrieval of geophysical parameters based upon remotely sensed data 
requires a sensor response model which relates the upwelling radiances that the sensor 
observes to the parameters to be retrieved. In the retrieval of precipitation water contents 
from satellite passive microwave observations, the sensor response model has two basic 
components. First, a description of the radiative transfer of microwaves through a 
precipitating atmosphere must be considered, because it is necessary to establish the 
physical relationship between precipitation water content and upwelling microwave 
brightness temperature. Also the spatial response of the satellite microwave sensor (or 
antenna pattern) must be included in the description of sensor response, since precipitation 
and the associated brightness temperature field can vary over a typical microwave sensor 
resolution footprint 

Olson (1989) utilized a simple cuboidal, radiative transfer cloud model to describe 
the upwelling brightness temperatures at the Scanning Multichannel Microwave Radiometer 
(SMMR) frequencies. Upwelling brightness temperatures were convolved by approximate 
radiometer antenna patterns to simulate the antenna temperatures measured by the SMMR. 
More recently, in applications to aircraft microwave data, Kummerow, et al. (1989) 
calculated upwelling microwave brightness temperatures for a collection of vertical 
precipitation profiles, allowing for variations in the total vertically-integrated precipitation 
by scaling the selected profile. In applications of this method to Special Sensor 
Microwave/Imager (SSM/I) data, vertically integrated precipitation water contents were 
assumed to vary lognormally in space within the raining fraction of the radiometer footprint 



(Kummerow, 1992). Sensor resolution effects were accomodated by preprocessing the 
SSM/I data using an antenna pattern matching technique (Robinson, et al., 1992). Xiang, 
et al. (1992) utilized three-dimensional cloud model simulations as the basis for upwelling 
brightness temperature calculations at the SSM/I frequencies. The retrieved profiles of 
cloud and precipitation profiles were those associated with modeled microwave brightness 
temperatures most consistent with SSM/I observations. An antenna pattern matching 
technique was utilized to account for the differing spatial resolution of the SSM/I channels; 
ref. Farrar and Smith (1992). 

In the present study, a "population" of convective cells, as well as stratiform clouds, 
are simulated using a computationally-efficient multi-cylinder cloud model. Ensembles of 
clouds selected at random from the population, distributed over a 25 km x 25 km model 
domain, serve as the basis for radiative transfer calculations of upwelling brightness 
temperatures at the SSM/I frequencies. Sensor spatial response is treated explicitly by 
convolving the upwelling brightness temperatures by domain-integrated SSM/I antenna 
patterns. The sensor response model is utilized in precipitation water content retrievals, 
which is the subject of Part II in this series. 

2. DEFINITIONS 

The sensor response model describes the physical relationship between the values of 
geophysical parameters P in the earth-atmosphere system, and the antenna temperatures 
TA which the sensor measures. A single antenna temperature measurement TA may 
respond to variations in several geophysical parameters at different locations; therefore the 
response model is written TA(P) m od to represent this multivariate dependency. 

It is assumed that the field of geophysical parameters P over the sensor swath can be 
represented by an array of values at discrete grid locations. The gridded domain for SSM/I 
retrievals is illustrated in Fig. 1. Rectangular grid "boxes" are defined on the grid, such 
that the center of each box coincides with the earth location of an all-channel antenna 


temperature measurement from the SSM/I. The boxes are oriented along the SSM/I A- 
scans (which contain the all-channel measurements) at a regular spacing of 25 km. Since 
the cross-scan separation of succeeding A- scan lines varies along the scan line, the cross- 
scan dimension of the grid boxes also varies, with a maximum dimension of 25 km in the 
center of the swath, and decreasing towards swath edge. In the discretization of the 
antenna patterns (Section 3) and brightness temperature fields (Section 4), the variation of 
the grid box dimension with scan position is taken into account. 

In the grid system depicted in Fig. 1 , an element of TA(P) mo d is defined by 

TA(P) p = A c T TB(P) p + A x T TB(P) p ' + SpTbb, (1) 

where A c and A x are the co- and cross-polarized antenna patterns of the radiometer, 
TB(P)p and TB(P)p' are the modeled brightness temperatures in the polarization p and 
orthogonal polarization p' with respect to the plane of polarization of the measurement, and 
5p is the fraction of the radiometer feedhom pattern not subtended by the antenna. Tbb is 
the cosmic background radiance (2.7 K). Each antenna temperature is modeled as the 
convolution of the upwelling brightness temperature field by the spatial response pattern of 
the sensor (see section 3 below). In this way, measurements at different channel 
frequencies or from different sensors which have different sampling / spatial resolution can 
be accommodated. 

In the following sections, the components of the sensor response model, Eq. (1), are 
described. 

3. SENSOR SPATIAL RESPONSE FUNCTIONS 

SSM/I measurements are diffraction -limited because the channel wavelengths (1.55 
cm to 0.35 cm) are not small in comparison to the antenna dimensions (61 cm x 66 cm); 
from Hollinger (1991). Aperture diffraction effects cause the breadth of the sensor antenna 


patterns to decrease with increasing frequency. The result is a large variation in the 
effective spatial resolution of the SSM/I. 

Measurements of the SSM/I antenna patterns and a detailed description of the SSM/I 
scan geometry were provided by Poe (1988). Based upon this information, the co- and 
cross-polarized SSM/I antenna patterns at the earth's surface are calculated for each SSM/I 
scan position and channel. The patterns are then area-integrated over each grid box in the 
swath grid and then normalized over all grid boxes to yield spatial response functions. 
Thus, the sensor spatial response functions A c and A x give the mean response, for each 
SSM/I measurement, to brightness temperatures upwelling from each grid box in the swath 
grid of Fig. 1. Representative normalized spatial response functions at 19.35, 37, and 
85.5 GHz are illustrated in Fig. 2. It may be noted from the figure that the sensor spatial 
response at 19.35 GHz is spread over several grid boxes in the swath grid, indicating 
relatively low spatial resolution at the grid scale (25 km). In contrast, the spatial response 
at 85.5 GHz is concentrated almost entirely in a single grid box, indicating high spatial 
resolution. The spatial response patterns are included in the sensor response model to 
accomodate the disparity in spatial resolution exhibited by the different channels of SSM/I. 

4. CLOUD ENSEMBLE/RADIATIVE MODEL 

The relationship between the geophysical parameters P and the brightness 
temperatures TBp and TBp' upwelling from an individual grid box is described using a 
cloud ensemble/radiative model. The model geometry is illustrated in Fig. 3. The 
horizontal dimensions of the domain are 25 km x 25 km, the nominal grid box dimension, 
with a vertical dimension of 20 km. Typically, within a region of these dimensions, an 
array of convective clouds as well as anvil-type or stratiform cloud may coexist. Both 
types of cloud are simulated in this study, with a brief summary given below. 


4. 1 Model Description 

Individual cloud simulations are performed using two versions of the multi-cylinder 
cloud model of Ridout (1991, 1993). Ridout's original model (ref. 1991, 1993) is utilized 
for convective cloud simulations. It is a time-dependent, non-hydrostatic cloud model 
which utilizes six concentric, quasi-cylindrical regions to represent the horizontal exchange 
of heat and moisture between the core region of the modeled convective cell and its 
immediate environment. The vertical momentum equation is represented explicity only in 
the innermost cylinder (core) of the cloud model, thereby reducing computational 
requirements substantially. The model includes a parameterization for the effects of vertical 
wind shear on the entrainment/detrainment of liquid water and other cloud variables. 

The treatment of microphysics follows the development by Lin, et al. (1983), with 
some modifications by Ferrier (1988). The model predicts the distributions of cloud water 
and cloud ice, rain, snow, and graupel/hail at levels separated by 300 m in the vertical, 
averaged over each cylindrical region. 

A model cloud is initially forced by a prescribed sub-cloud vertical velocity 
distribution over a prescribed period (900 s in the current study), after which time the 
forcing is removed and the cloud continues to develop according to model dynamics. This 
method of initiation is used to simulate gust-front forcing of convective clouds. Cloud 
development is also controlled by a prescribed minimum updraft radius which is the radial 
dimension of the updraft (third cylinder) at midlevels in the troposphere. The radii of the 
other cylinders are all scaled according to the prescribed updraft radius. 

A second version of the Ridout (1992) model is developed to simulate stratiform 
anvil-type precipitation distributions in this study. In this version of the model, the five 
surrounding cylinders are eliminated, and it is assumed that the distributions of all cloud 
variables both inside and outside the remaining cylinder are identical. In other words, a 
strictly one-dimensional simulation is performed. In addition, the vertical velocity 
distribution is prescribed to be a double-parabolic profile, with a downdraft below the 


freezing level and an updraft aloft. This type of vertical velocity profile is generally 
consistent with observations documented by Leery and Houze (1979), Gamache and 
Houze (1982), Marks and Houze (1987), and Rutledge and Houze (1987), who observed 
mesoscale updrafts and downdrafts in several tropical and midlatitude storms. The 
Lin/Ferrier cloud microphysical parameters are altered to reflect the variability of particle 
size distributions noted by Houze, et al. (1979) and Stewart, et al. (1984). Stratiform 
cloud simulations are performed by integrating the model forward in time until equilibrium 
distributions of all cloud variables are obtained, usually after 10 hours of simulation time. 

4.2 Cloud Environment 

The cloud environment for both convective and stratiform cloud simulations is taken 
from the GATE Day 261 rawinsonde profile; ref. Ferrier and Houze (1989). Profiles of 
pressure, temperature, water vapor mixing ratio, and horizontal wind components are 
interpolated to levels at 150 m altitude spacing in the Ridout model. 

4.3 Generation of Ensembles - Convective Clouds 

By varying the minimum updraft radius and the basal updraft vertical velocity in the 
Ridout model, a variety of clouds ranging from shallow to deep convection are simulated. 
Combinations of eight cloud radii between 0.5 and 5 km and four basal updrafts between 2 
and 8 m/s are used to generate 32 different cloud simulations. For each cloud simulation an 
effective cloud lifetime is computed. For the purpose of this study the lifetime of a 
simulated cloud is defined as the period over which the average reflectivity of the cloud 
core volume (a cylindrical region 2.43 km in radius and 9.49 km in height) exceeds 30 
dBZ. Effectively, these are the cloud radar echo lifetime criteria that were used by L6pez 
(1977) in his study of tropical clouds. The core volume dimensions are designed to 
simulate the radar volume of the airborne APS-20 radar utilized in Ldpez' study. 

Based upon the radar echo lifetime statistics compiled by L6pez (1977), a 
"population" of simulated clouds is created. In Ltipez' study, tropical clouds were found to 
have echo lifetimes that were lognormally distributed. A plot of this observed lifetime 


distribution is reproduced in Fig. 4. A simulated cloud with a given lifetime is assumed to 
occur with the same relative frequency as a cloud with the same observed lifetime in 
Ldpez' radar study. 

A random number generator is utilized to select simulated clouds from the 
population and create an ensemble within the 25 km x 25 km model domain (see Fig. 3). 
Once a cloud is selected from the population, a cloud lifestage is also randomly selected. 
The cloud from that stage of the simulation is then positioned at random within the model 
domain. Another cloud is then selected at random from the assumed population, and the 
procedure is repeated. The only restriction on cloud placement is that the updrafts of no 
two simulated clouds may overlap. By varying the number of clouds selected for the 
ensemble, a wide range of cloud configurations with different water and ice contents can be 
simulated. In this study, 350 different cloud ensembles containing between 8 and 56 
clouds are created using the procedure outlined above. One such ensemble based upon 48 
simulated clouds is depicted in Fig. 5. In the figure, image intensity is proportional to the 
slant-path integrated water contents at the SSM/I viewing incidence angle (53.1 degrees 
from zenith). 

4.4 Generation of Ensembles - Stratiform Clouds 

Different stratiform clouds are generated by varying the amplitude of the prescribed 
updraft and downdraft distributions. Following the diagnostics studies of Gamache and 
Houze (1982), in which the magnitudes of updraft and downdraft velocities were estimated 
for tropical anvils, a range of representative updraft and downdraft amplitudes are selected. 
Updraft amplitudes are varied from 0.05 to 0.50 m/s, while downdraft amplitudes are 
varied from 0.00 to -0.50 m/s, to create a total of 42 stratiform simulations. 

4 . 5 Radiance Calculations 

Both convective and stratiform cloud ensembles are located within the model domain 
depicted in Fig. 3. The 25 km x 25 km x 20 km model domain is subdivided into a grid 
with a horizontal grid spacing of 0.500 km, and a vertical grid spacing of 0.375 km. 


Pressure, temperature, water vapor density, and the equivalent water contents of 
precipitating and non-precipitating particles from the simulated clouds are interpolated to 
grid locations within the domain. In addition, cloud latent heating rates are evaluated at 
grid locations. The average surface rainfall rate over the domain is also computed. 

The cloud and environmental parameters at each grid location are used to calculate the 
associated microwave radiative properties (i.e., the extinction coefficients, single-scatter 
albedoes, and asymmetry factors). Molecular oxygen and water vapor extinction 
coefficients are computed using the Liebe (1985) model. Microwave radiative parameters 
for precipitating liquid and ice hydrometeors are computed using the formulae of 
Kummerow and Weinman (1988). Rayleigh theory is employed to calculate the extinction 
coefficients of nonprecipitating cloud particles. At the base of the grid a surface skin 
temperature and microwave emissivities characteristic of both land and ocean surfaces are 
assigned. The downwelling cosmic background brightness temperature (2.7 K) is 
specified at the top of the model domain. 

A forward radiative transfer calculation is performed to determine the distribution of 
microwave brightness temperatures upwelling from each cloud ensemble domain. A 
radiative transfer solution based upon Eddington's second approximation is utilized to 
accomodate the multiple scattering effects of precipitating hydrometeors at microwave 
frequencies. Brightness temperatures are computed at the SSM/I frequencies (19.35, 
22.235, 37, and 85.5 GHz) in both vertical and horizontal polarizations at an incidence 
angle of 53. 1 degrees. Since radiance paths at oblique incidence can exit the sides of the 
model domain, it is assumed that the domain, and all cloud and environmental properties, 
are periodic in the horizontal. 

Fields of brightness temperatures upwelling from the convective cloud ensemble of 
Fig. 5 are depicted in Fig. 6 (over a low-emissivity background characteristic of the ocean), 
and in Fig. 7 (over a high-emissivity backgound characteristic of a land surface). 
Comparing the figures, it may be noted that the ocean background provides a good 



radiative contrast to regions of greater emission from liquid precipitation, while the higher- 
emissivity land background provides less contrast Over either background, the low 
brightness temperatures at 85.5 GHz associated with scattering by precipitating ice particles 
are evident. The basic differences between liquid and ice precipitation signatures provide 
the motivation for the radiative parameterization described in the next section. 

Ensemble- average, horizontally-polarized brightness temperatures are plotted against 
the corresponding ensemble-average slant path-integrated precipitating liquid water contents 
for the four SSM/I frequencies in Fig. 8. The brightness temperature calculations in these 
plots utilize an emissivity of 0.3, characteristic of an ocean surface. Both convective and 
stratiform ensemble calculations are included in the figure. Plotted numbers in the figure 
are the ensemble-average slant-path integrated ice water contents, thresholded at integral 
values. 

It may be noted from the figure that there is a significant spread of brightness 
temperatures corresponding to a given liquid precipitation path, due to variations in ice 
water content, precipitation area, non-precipitating cloud amount, and cloud geometry. The 
scattering effect of ice precipitation is most obvious at 37 and 85.5 GHz, leading to 
somewhat lower brightness temperatures for a given precipitating liquid path. At 85.5 
GHz, the stratiform ensembles yield brightness temperatures which are significantly lower 
than the convective ensemble brightness temperatures, due to the relatively high amounts of 
precipitating ice (relative to precipitating liquid) in the stratiform clouds. The focus of the 
next section will be to parameterize the ensemble-average upwelling brightness 
temperatures as functions of the precipitating liquid and ice paths. 

5 . RADIATIVE PARAMETERIZATION OF CLOUD ENSEMBLE/ 

RADIATIVE CALCULATIONS 

The cloud ensemble/radiative calculations provide a fairly detailed and accurate 
description of the upwelling brightness temperature field for a given cloud ensemble 
configuration. However, the sensor response model (Eq. 1) requires a functional 


relationship between the grid box average brightness temperatures (TB p and TB p ') and 
properties of the cloud field and its environment (geophysical parameters P). This 
functional relationship, or parameterization, must include a meaningful description of the 
precipitation field and produce an accurate representation of the detailed cloud 
ensemble/radiative calculations. However, since the sensor response model will ultimately 
be utilized in retrievals, the number of geophyical parameters P must be consistent with the 
sampling of the SSM/I. The SSM/I provides only 13 brightness temperture measurements 
per grid box in the swath grid of Fig. 1., and it is likely that some of the information in 
these measurements is redundant. Thus the number of parameters must be kept small in 
order to avoid ambiguous retrievals; i.e. the number of unknown parameters P in the 
retrieval must be less than the number of independent measurements. Our strategy is to 
select parameters which describe the mean cloud and environmental characteristics of each 
grid box, and then parameterize smaller scale variability. 

Figures 5 through 7 illustrate the basic relationships between precipitating liquid and 
ice contents and upwelling microwave brightness temperatures. There is a strong 
correlation between the path-integrated precipitating water content and the upwelling 
brightness temperature at 37 GHz. The same correlation is observed in the 19.35 and 
22.235 GHz simulations (not shown). At 85.5 GHz, the path-integrated precipitating ice 
content is strongly anticorrelated with upwelling brightness temperature. 

These relationships suggest a simple parameterization for the average upwelling 
brightness temperature from a given swath grid box. The full three-dimensional radiative 
calculation is replaced by four one-dimensional calculations. The first calculation is 
performed for the brightness temperature upwelling from the cloud-free environment of the 
grid box. The second assumes that only precipitating liquid is present in the model 
atmosphere; the third only precipitating ice, and the fourth for a model atmosphere 
containing both precipitating liquid and ice. Vertical profiles of precipitating liquid and ice 
in the one-dimensional radiative calculations are represented by analytical curve fits to the 



horizontally-averaged profiles of both the convective and stratiform cloud ensembles. Both 
precipitating liquid and ice equivalent water contents are given by 


| 0, z > zmax 


PLWC 


A • sin (k (z - zmin)/(zmax - zmin)) (2) 

• exp (-{Bl + B2/PLWP) • (z - zmin)/(zmax - zmin)) , 


zrrnn < z < zmax 


0, z < 0 

Here, PLWC is the precipitating liquid or ice equivalent water content, PLWP is the slant 
path-integrated precipitating liquid or ice depth, zmin and zmax are the lower and upper 
altitude bounds of the water content profile, and Bl and B2 are fitting constants. Given the 
slant path-integrated precipitating liquid water or ice depth, the constant A is set such that 
the integrated analytical profile yields the same depth. The constants zmin, zmax, Bl, and 
B2 are all optimized to yield the minimum mean-square difference between the water 
content profiles generated by the cloud ensembles and those calculated using Eq. (2). 
Values of the adjustable parameters and statistics of the analytical profile approximation are 
listed in Table 1. 

The mean profiles and analytical curve fits for convective liquid and ice precipitation 
are shown in Figs. 9 and 10, respectively. Similar curve fits are obtained for the stratiform 
cloud ensembles. The mean error standard deviation in the estimated water contents from 
the analytical profiles is about 0.03 g/m^. 

The four one-dimensional calculations are combined in proportion to the fractional 
coverage of liquid and ice precipitation that are observed along the oblique viewing angle of 
the SSM/I. Although it is assumed that precipitating liquid and ice occupy specified 
fractions, no assumption is made concerning the relative positions of the liquid and ice 


fractions with respect to one another inside the grid box. If fj and fi are the fractions of 
precipitating liquid and ice, respectively, then the parameterized upwelling brightness 
temperature at a given frequency and polarization is given by 

TB(P) = (1 - fi) (1 - fi) TB C + fl (1 - fi) TBi 

+ (1 - fl) fi TBi + fl fi XBii . (3) 

Here, TBc, TB], TBi, and TBu are the brightness temperatures upwelling from a cloud-free 
model atmosphere, an atmosphere containing only precipitating liquid, an atmosphere 
containing only preciptating ice, and an atmosphere containing both precipitating liquid and 
ice, respectively. 

Although Eq. (3) is a fairly good approximation to the brightness temperatures 
upwelling from the stratiform cloud ensembles, the small-scale variability of precipitation 
water paths in the convective ensembles leads to discrepancies which cannot be 
parameterized so easily. To account for the variability of water paths, the presence of 
nonprecipitating clouds, and other factors, the liquid and ice precipitation water contents in 
the analytical profiles are adjusted to obtain better agreement with the detailed brightness 
temperature calculations. This strategy follows the modeling efforts of Short and North 
(1990), and Hinton, et al. (1992), who noted that the brightness temperature response to 
the average precipitation rate is less if the rain rates are gamma or lognormally distributed, 
rather than uniformly distributed. Thus, the effective liquid and ice precipitation water 
contents are expressed as power-law functions of the original water contents (from the 
analytical profiles), to accomodate the difference in response. 


PLWC' = a PLWC P 


( 4 ) 


Here PLWC is the original liquid or ice water content. The power-law constants a and (3 
are adjusted to minimize the difference between the parameterized brightness temperatures, 
Eq. (3), and those computed from the cloud ensemble/radiative calculations. Separate 
power-law constants are computed for the convective and stratiform ensembles. Best-fit 
power law constants and parameterization errors are summarized in Table 2. Also 
scatterplots of the parameterized brightness temperatures vs. the cloud ensemble brightness 
temperatures at the four SSM/I frequencies are presented in Fig. 1 1. 

It may be noted from Table 2 that the power-law constants used to modify the 
convective precipitating liquid water contents at 19.35, 22.235, and 37 GHz are all 
significantly less than 1, indicating a substantial reduction in the effective water content 
utilized in the radiative parameterization. This result is consistent with the study of Hinton, 
et al. (1992), who noted a reduced response of up welling brightness temperature to mean 
rainfall rate if the rainfall rates within a given region are gamma-distributed. The power- 
law constants used to modify the stratiform precipitating liquid water contents are closer to 
1, which is reasonable since the parameterization appoximates the precipitation field by a 
horizontally-uniform field, and the stratiform cloud ensembles have horizontally-uniform 
structure. 

The power-law constants for convective ice precipitation defy a simple physical 
interpretation; however, the power-law constants used to modify the stratiform precipitating 
ice equivalent water contents are close to 1 . Again, the horizontally-uniform clouds in the 
radiative parameterization are a good approximation to the true structure of the stratiform 
cloud ensembles, so only a slight modification of the water contents is required. 

The statistics in Table 2 and the scatterplots in Fig. 1 1 indicate parameterization errors 
which are typically small, but which increase with increasing frequency. It should be noted 
that the absolute accuracy of SSM/I antenna temperature measurements is on the order of 
± 3 K (Hollinger, 1991). The most serious parameterization errors are attributed to the 
strongly convective ensembles which produce brightness temperatures between 180 K and 


200 K at 85.5 GHz (see Fig. lid). Scattering by precipitation-sized ice particles is 
strongest for these convective clouds at 85.5 GHz, and the horizontally-homogeneous 
cloud approximation utilized in the parameterization is least suitable under these conditions. 

In the retrieval scenario it will not be known a priori whether a given brightness 
temperature is upwelling from a convective or stratiform cloud field. Therefore, a simple 
discriminant based upon the relative proportions of liquid and ice precipitation is devised. 
The distribution of ice vesus liquid preciptation paths of the convective and stratiform 
ensembles is shown in Figure 12. It may be noted that the ratio of ice to liquid paths is 
much smaller in the convective ensembles than in the stratiform ensembles. In the 
convective ensembles, the ratio of ice precipitation path to the total path, designated rj ce , is 

always less than 0.5, while in the stratiform ensembles it typically exceeds 0.8. In this 
study it is assumed that if rjce is less than 0.5, the convective parameterization for 
brightness temperature is used. If ri ce is greater than 0.8, then the stratiform 
parameterization is employed. For ratios between 0.5 and 0.8, it is assumed that the 
precipitation water paths have both convective and stratiform components, the proportion 
of each being determined by linear interpolation. 

6. CONCLUDING REMARKS 

The cloud ensemble/radiative parameterization described herein provides a 
computationally-efficient means of simulating microwave antenna temperatures at the 
SSM/I frequencies. The parameterization is therefore useful in physical retrieval 
applications where many thousands of antenna temperature computations must be 
performed over a single solution domain (see Pan II of this series). The partitioning of 
cloud fields into separate convective and stratiform components is an obvious 
simplification, since generally speaking there are dynamical and thermodynamic 
interactions between the convective and stratiform regions of an organized precipitating 
storm; ref. Rutledge and Houze (1987). The cloud ensemble/radiative parameterization 



could be generalized by simulating cloud fields in a fully 3-dimensional cloud model and 
parameterizing the bulk precipitation structure and radiance fields as in the current study. 

Other data from the SSM/T, SSM/T-2, and infrared sensors could be included in 
physical retrievals by specifying the appropriate sensor spatial response characterisitics and 
channel frequencies in the cloud ensemble/radiative parameterization. The addition of these 
sensor data will be the subject of future study. 
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Table 1. Parameters of the best-fit analytical precipitating water content profiles, based 
upon the convective and stratiform cloud ensemble simulations. Adjustable parameters 
zmin and zmax are in units of kilometers, B1 is dimensionless, and B2 has units of m 2 /kg. 
Parameters associated with the liquid and ice precipitation profiles are listed separately. 

The analytical convective profiles have a combined bias of -2.6 x 10*5 g/m^, and an error 
standard deviation of 0.037 g/m^. The analytical stratiform profiles have a combined bias 
of 1.1 x 10"4 g/m^, and and error standard deviation of 0.028 g/m^. 



zmin 

zmax 

B1 

B2 

Convective 

liquid: 

-0.84 

6.19 

1.81 

1.67 

Convective 

ice: 

3.75 

15.20 

5.25 

3.22 

Stratiform 

liquid: 

-0.77 

5.07 

0.99 

-0.29 

Stratiform 

ice: 

3.13 

10.70 

3.21 

-0.37 


™e2 Best-fit power-law constants in the cloud ensembleAadiative parameterization at 
the SSM/I frequencies. Constants are listed separately for convective and stratiform cloud 
ensembles. Also listed are the bias and error standard deviation of the parameterized 
brightness temperatures with respect to the brightness temperatures generated from the 
cloud ensembles. 
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Fig. 1. Schematic of the SSM/I swath grid. For clarity, grid box sizes have been 
magnified by a factor of three. 
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Fig. 2. Co-polarized SSM/I response functions at (a) 19.35 GHz, (b) 37 GHz, and (c) 
85.5 GHz in the vertical polarization channels. 









Fig. 4. Radar echo lifetime cumulative distribution from Ldpez (1977). The cloud radar 
echo lifetime tlife, in minutes, is plotted against the cumulative fraction of clouds F with 
lifetimes less than tlife. 
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Fig. 5. Imagery of (a) path integrated precipitating liquid water and (b) path integrated 
precipitating ice from the cloud ensemble/radiative model. The precipitation path 
distributions are generated from an ensemble of 48 convective cloud simulations. The area 
represented by each panel is 25 km x 25 km. 



Fig. 6. Computed brightness temperature field at (a) 37 GHz, and (b) 85.5 GHz, 
upwelling from the convective cloud ensemble shown in Fig. 5. A surface emissivity of 
0.30 is assumed. The area represented by each panel is approximately 25 km x 25 km. 



Fig. 7. Computed brightness temperature field at (a) 37 GHz, and (b) 85.5 GHz, 
upwelling from the convective cloud ensemble shown in Fig. 5. A surface emissivity of 
0.95 is assumed. The area represented by each panel is approximately 25 km x 25 km. 
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Fig. 8. Plot of cloud ensemble/radiative model brightness temperatures vs. ensemble- 
average precipitating liquid water paths at (a) 19.35 GHz and (b) 22.235 GHz. Plotted 
numbers are ensemble-average precipitating ice water paths thresholded at integral values. 







•5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

LIQUID PRECIPITATION PATH [KG/H«2] 


DU CT 0 0 u 

°Vo°° l 
S° * ' \ 


\, 0 .5 l.o 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

\J LIQUID PRECIPITATION PATH IKG/M..21 

Fig. 8. (corn.) Plot of cloud ensemble/radiative model brightness temperatures vs. 
ensemble-average precipitating liquid water paths at (c) 37 GHz and (d) 85.5 GHz. Plotted 
numbers are ensemble-average precipitating ice water paths thresholded at integral values. 
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Fig. 9. Mean profiles of precipitating liquid water from (a) the convective cloud ensemble 
simulations, and (b) the analytical curve fits described in the text. Profiles are shown for 
precipitating liquid paths of 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 kg/m 2 . 
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Fig. 10. Mean profiles of precipitating ice from (a) the convective cloud ensemble 
simulations, and (b) the analytical curve fits described in the text. Profiles are shown for 
precipitating ice paths of 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 kg/m 2 . 
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Fig. 11. (cont.) Plot of parameterized brightness temperatures vs. cloud ensemble/radiative 
model brightness temperatures at (c) 37 GHz and (d) 85.5 GHz. For reference, the solid 
line in each panel represents a perfect correlation. 
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Fig. 12. Scatterplot of ice precipitation path versus liquid precipitation path from both the 
convective cloud ensemble (C) and stratiform cloud ensemble (S) simulations. Ice to total 
precipitation path ratios of 0.5 and 0.8 are indicated by solid lines. 
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1. INTRODUCTION 

The retrieval of precipitation distributions by passive microwave techniques has been 
an area of study in satellite remote sensing for well over a decade. The recent initiative to 
estimate both the horizontal and vertical distributions of precipitating liquid water in the 
tropics, i.e. the Tropical Rainfall Measuring Mission or TRMM, has further stimulated 
interest in this area; see Simpson, et al. (1988). The proposed TRMM satellite will carry a 
suite of sensors including passive and active microwave instruments, as well as visible and 
infrared radiometers. In addition to TRMM, a combination of microwave and visible/IR 
sensors has also been proposed for the Earth Observing System (EOS) platforms for the 
purpose of retrieving rainfall rates and other geophysical parameters; ref. Schoeberl, et al. 
(1992). Methods for determining precipitation distributions based upon data from several 
instruments having different spectral and resolution characteristics will therefore be 
required. 

Currently multispectral microwave observations from the Special Sensor 
Microwave/Imager (SSM/I), borne by the DMSP-F8, F10, and Fll, have provided 
researchers with the means for retrieving rainfall rates; see Hollinger (1991), Kummerow, 
et al. (1989), Petty and Katsaros (1990), Xiang, et al. (1992). The SSM/I is a 
multichannel, dual-polarization, passive microwave radiometer with channels at 19.35, 
22.235, 37, and 85.5 GHz. The diffraction limitation of the SSM/I antenna causes the 
spatial resolution of measurements to increase with frequency, such that the minimum 
footprint dimension decreases from 43 km at 19.35 GHz to 13 km at 85.5 GHz. The 
DMSP-F8, F10, and Fll also carry the Special Sensor Microwave/Temperature (SSM/T) 
sounder with channels at 50.5, 53.2, 54.35, 54.9, 58.4, 58.825, and 59.4 GHz, and a 
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footprint dimension of approximately 180 km at nadir view. In addition the DMSP FI 1 
carries the Special Sensor Microwave/Temperature -2 (SSWT-2) sounder, with channels at 
91.665, 150, and 183.3 ± 1, ± 3, ± 7 GHz, and a minimum footprint dimension of 48 km 
at nadir view. The SSM/T and SSM/T-2 can provide temperature and humidity information 
not available from the SSM/I alone. 

In the present study, a physical retrieval method for estimating precipitating water 
distributions and other geophysical parameters based upon measurements from the DMSP- 
F8 SSM/I is developed. Three unique features of the retrieval method are (1) sensor 
antenna functions are explicitly included to accommodate varying channel resolution, (2) 
precipitation-brightness temperature relationships are quantified using the cloud 
ensemble/radiative parameterization described in Part I, and (3) spatial constraints are 
imposed for certain background parameters, such as humidity, which vary more slowly in 
the horizontal than the cloud and precipitation water contents. The general framework of 
the method will facilitate the incorporation of measurements from the SSM/T, SSM/T-2 and 
geostationary infrared measurements, as well as information from conventional sources 
(e.g. radiosondes) or numerical forecast model fields. 

2. METHOD 

The basis of the retrieval method is the minimum variance solution as described in 
Lorenc (1986). An error functional which expresses the deviation of the observed satellite 
antenna temperatures from model-derived values, plus an additional constraint which 
represents the deviation of retrieved geophysical parameters from a priori estimates, is 
derived. 

E = (TA 0 bs * TA(P) m od) T ^ 1 (TA 0 bs * TA(P) mo d) 

+ (P - Papriori)^ (P-Papriori) 



( 1 ) 


Here, TA 0 bs is a vector of antenna temperature measurements, TA(P) mo d are antenna 
temperature estimates based on the sensor response model (see Part I of this series), P is a 
vector of geophysical parameters to be retrieved, and Papriori are estimates of the 
parameters based upon a priori knowledge. E and A are covariance matrices representing 
observational plus model errors (L), and errors in the a priori estimates (A). In the present 
study all off-diagonal elements of these matrices are neglected. 

The retrieval of precipitating liquid water contents and other geophysical parameters is 
accomplished by iteratively perturbing the geophysical parameters until the error functional 
Eq. (2) is minimized. Unlike Lorenc (1986), the error functional is minimized using the 
memoryless, quasi-Newton method described by Shanno (1978). The quasi-Newton 
approach is computationally efficient when the retrieval problem requires the simultaneous 
optimization of a large number of unknown geophysical parameters. In this situation, 
matrix inversion methods are impractical. 

3. APPLICATION OF THE RETRIEVAL METHOD 

TO SSM/I DATA 

3. 1 Retrieval Parameters and Constraints 

In Part I of this series, a sensor response model was developed for SSM/I to relate 
measured microwave antenna temperatures to variations in geophysical parameters in the 
earth/atmosphere system. The sensor response model is designated T A (P)mod in Eq- (1). 
The explicit dependence of antenna temperature on slant path-integrated liquid and ice- 
phase precipitation, fractional coverage of liquid and ice precipitation within a 25 km x 25 
km grid box, atmospheric temperature and water vapor, and surface temperature and 
emissivity were discussed in Part I. Although each of these geophysical parameters can 
have a significant effect on SSM/I antenna temperatures, only a limited number of 
unknown parameters can be allowed in the minimization of the error functional Eq. (1). 
Otherwise, ambiguous solutions can result. 



Of primary interest are the distributions of precipitating liquid and ice over the swath 
grid, and therefore the slant path-integrated preciptitating liquid and ice, and the associated 
grid-box fractions are allowed to vary in the retrievals. Although variations in surface skin 
temperature and emissivity can affect microwave surface emission, climatological variations 
of surface skin temperature are relatively small. Surface skin temperature over the ocean 
varies by less than ± 5 K in the tropics, up to about ± 10 K at midlatitudes (Pickard and 
Emery, 1982). Any error in the retrieved value of skin temperature can easily be 
compensated by a change in the surface emissivity to produce nearly the same surface 
emission. Therefore, in this study the value of surface skin temperature is prescribed at a 
climatological value representative of the tropics (300.15 K). A somewhat better estimate 
could be obtained from infrared retrieval composites, but a climatological value is found to 
be sufficiently accurate for the purpose of demonstrating the retrieval method. Surface 
emissivities are parameterized according to the model of Grody (1988). The emissivity in 
each polarization has a frequency dependence given by 

e = a + b * In (v) (2) 

where v is the channel frequencey and a and b are free parameters. Separate emissivity 
parameters describing the emissivities of land and ocean surfaces are allowed to vary in the 
retrieval, but it is assumed that these parameters are the same for every grid box of the 
swath grid. To account for the mixture of emissivities in grid boxes which contain both 
land and water surfaces, an additional land fraction parameter, f] a nd» is also allowed to 
vary in the retrieval. The effective emissivity is given by 


eeffective = fland * e land + (1 * fland) * e ocean 


(3) 


In the retrieval method, profiles of atmospheric temperature and water vapor mixing 
ratio are represented by a mean profile and a set of profile eigenfunctions. The mean 
profiles and eigenfunctions are derived from a global radisonde data base of 1200 profiles 
developed by Smith, et al. (1974). 

A link between the atmospheric profiles and cloud ensemble/radiative 
parameterization is established as follows. A profile of temperature and humidity (water 
vapor mixing ratio) is selected to serve as the environment for the cloud model simulations 
and ensemble generation described in Part I. The relationship between the generated 
ensembles and the computed upwelling microwave brightness temperatures is 
parameterized, and a set of fitting constants is calculated for that particular environment (see 
also Part I). Additional profiles of atmospheric temperature and humidity can be selected 
to generate cloud ensemble/radiative parameterizations appropriate for other climatologies. 

In the retrieval, an initial guess profile is compared to all of the profiles utilized to 
generate cloud ensemble/radiative parameterizations. The parameterization associated with 
the profile which best matches the initial guess profile (in the mean-square sense) is 
selected for the retrieval of precipitation amounts in a given grid box. The guess profiles 
are allowed to deviate from the parameterization profiles by allowing the weighting of the 
first eigenfunctions of temperature and mixing ratio (in an eigenfunction decomposition of 
the guess profiles) to vary in the retrieval. 

In the present study, the GATE Day 261 profile (ref. Ferrier and Houze, 1989) is 
selected to yield a cloud ensemble/radiative parameterization representative of tropical 
conditions. This single profile is selected for the demonstration of the precipitation retrieval 
method, since the applications herein include SSM/I observations of a tropical squall line 
and a developing tropical storm. However, cloud ensemble/radiative parameterizations 
corresponding to other temperature and humidity environments could be added to represent 
midlatitude conditions or other climatologies. 



Physical retrieval tests in which the temperature and mixing ratio eigenfunction 
weights were allowed to vary in each grid box resulted in a lack of constraint on solutions 
and ambiguous retrievals. Ambiguous solutions were due to the fact that temperature 
variations and mixing ratio variations over land produce only a very small signal in 
upwelling microwave brightness temperatures at the SSM/I frequencies, and these signals 
are generally indistinguishable from variations in background emissivity or precipitation 
amount. 

The ambiguities are removed by prescibing the first eigenfunction weight of 
temperature over the entire swath grid, and allowing the first eigenfunction weight of water 
vapor mixing ratio to vary only over oceanic regions. Over land, the first eigenfunction 
weight of mixing ratio is again prescribed. 

Regarding the retrieval of precipitation parameters, preliminary testing also revealed 
an ambiguity between slant path-integrated precipitating ice and the fraction of precipitating 
ice within each grid box. The ambiguity in these retrieved quantities is removed by setting 
the fraction of precipitating ice equal to the fraction of precipitating liquid. This constraint 
is likely to introduce some error into retrieved precipitation amounts, since precipitating ice, 
which is associated with stronger convective updrafts and organized mesoscale motions 
(see Part I), generally has an areal extent which is different from the areal extent of 
precipitating liquid. The ice content/fraction ambiguity might be removed with the addition 
of higher- frequency ice sensitive channels on future microwave sensors. 

A summary of parameters P in the SSM/I physical retrieval and the associated 
constraints on the parameters is summarized in Table 1. 

In the initial demonstration of the retrieval method which is documented here, only 
the cloud ensemble/radiative parameterization based upon the GATE Day 261 profile is 
considered. This parameterization is appropriate since only retrievals of precipitation in the 
tropics are considered here. Precipitation retrieval applications in other climatic 



environments, and the selection of alternative cloud ensemble/radiative parameterizations 
appropriate to those environments, will be the subject of a future study. 

Estimates of observational and cloud ensemble/radiative parameterization errors are 
specified in the diagonal elements of E of Eq. (1). The observational error standard 
deviations of each SSM/I measurement are taken to be 3 K as suggested by Hollinger, et al. 
(1991) for absolute calibration errors. Parameterization error standard deviations are 
assumed to be 6 K in the 85.5 GHz channels and 3.6 K in the remaining channels. These 
figures reflect the greater residual errors in the 85.5 GHz parameterization fits to the cloud 
ensemble/radiative calculations of brightness temperatures (see Part I). An additional 
weighting factor of 0.25 is applied to the 85.5 GHz terms in Eq. (1) to compensate for the 
SSM/I 85.5 GHz sampling, which is four times as great as the sampling in the other 
channels. 

The minimization of E in Eq. (1) requires the evaluation of finite difference 
derivatives of brightness temperatures with respect to the parameters P. Efficient 
convergence to a minimum requires that these derivatives be scaled to compensate for 
variations in sensitivity of brightness temperature to a particular parameter. Appropriate 
scaling factors for each parameter are selected based upon extensive retrieval tests. 

3.2 Darwin Case Study 

The physical retrieval method, Eq. (1), is applied to SSM/I observations of a tropical 
squall line at 0906 UTC on February 9, 1988 in the vicinity of Darwin, Australia. SSM/I 
observation of the squall line at 37 GHz and 85.5 GHz in the horizontal polarization are 
shown in Figs. 1 a and 1 b, respectively. In the panels of the figure, the northern coast of 
Australia and the outlines of Bathurst and Melville Islands are indicated by white lines, and 
the grayscale of brightness temperature values in K is displayed at the bottom. Each panel 
has dimensions of approximately 240 km x 240 km. The convective leading edge of the 
squall line parallels the coast of Australia, and is indictated by a line of relatively low 
brightness temperatures in the 85.5 GHz channel. These depressed brightness 


temperatures are induced by the scattering of microwave radiances by precipitation-sized ice 
particles which are held aloft by relatively strong updrafts within the convective zone. A 
sequence of radar imagery from the Darwin/TOGA radar bracketing the SSM/I overpass 
time indicates a general southeast to northwest motion of the squall line (Steiner, 1992). A 
broad, trailing stratiform precipitation region to the south and east of the convective leading 
edge is seen as a relatively uniform depression in the 37 and 85.5 GHz brightness 
temperatures extending beyond the lower edge of each panel. The brightness temperatures 
in the trailing stratiform region tend to be higher than in the convective zone, which 
suggests that the concentrations of scattering ice particles are not as high in the stratiform 
zone. Regions of smaller-scale convective precipitation are seen as scattering depressions 
in the 85.5 GHz imagery to the northwest of the squall line, between Bathurst and Melville 
Islands, and also to the west of the convective line. 

The SSM/I retrieval method is applied to a 64 x 64 section of the swath grid centered 
on the Darwin/TOGA radar site, located at the center of the panels in Fig. 1 . All seven of 
the SSM/I channels are included in the retrieval. After 21 iterations of the quasi-Newton 
procedure the error functional E in Eq. (1) is effectively minimized without any additional 
constraints (second term equal to zero) to produce the retrievals of average slant path- 
integrated precipitating liquid and precipitating ice depicted in Figs. 2a and 2b, respectively. 
The image intensity in these panels is proportional to the equivalent column-integrated 
liquid water or ice from the retrieval. A grayscale relating depth-integrated water in kg/m^ 
to image intensity is shown at the bottom of each panel. 

From Fig. 2, it may be noted that the highest concentrations of integrated liquid 
precipitation occur along the convective leading edge of the squall line, with lower values 
spread almost uniformly over the trailing stratiform region. Peak values are close to 4 
kg/m^ in the convective leading edge, with typical values between 1 and 2 kg/m^ in the 
trailing stratiform region. Regions of liquid precipitation associated with the convective 
zones west and northwest of Darwin are also seen in the retrieval. Retrievals of the 



column-integrated precipitating ice in Fig. 2b exhibit relatively high values in both the 
convective leading edge and stratiform zones, on the order of 1 to 1.5 kg/m^. The 
retrievals in the convective and stratiform zones have roughly the same magnitude, in 
contrast to the retrievals of liquid precipitation, which exhibit greater values along the 
convective leading edge. Also, the maximum ice amounts in the stratiform zone tend to be 
concentrated to the southeast of Darwin. Precipitating ice amounts are significantly less in 
the convective regions northwest and west of Darwin, with maximum values of about 0.8 

kg/m^. 

A nearly concident radar volume scan from the 5-cm Darwin/TOGA radar starting at 
0850 UTC on February 9, 1988 is processed to obtain comparative estimates of integrated 
liquid and ice precipitation amounts. The radar reflectivities from the volume scan are first 
interpolated to a cartesian grid with overall dimensions 240 km x 240 km x 18 km. The 
grid spacings are 2 km in the horizontal and 1.5 km in the vertical. The radar data in 
elevation/azimuth/range format are interpolated to grid locations using the method of Mohr 
and Vaughan (1979). 

Within the convective leading edge of the squall line, the fractions of liquid and ice 
precipitation associated with a given radar reflectivity are estimated using the temperature- 
dependent partitioning scheme of Moss and Johnson (1992). Within the trailing stratiform 
region, it is assumed that the radar reflectivity is due entirely to ice-phase hydrometeors 
above the freezing level, and below the 6 C level the radar reflectivity is due entirely to 
liquid hydrometeors. A linear transition from ice to liquid between 0 C and 6 C is 
assumed. This stratiform precipitation partitioning scheme is based upon model 
simulations of anvil clouds using a one-dimensional version of the Ridout (1991, 1993) 
model. The vertical profile of environmental temperatures is obtained from the Darwin 
sounding at 2200 UTC on February 9, 1988. 

Slant path-integrated liquid and ice-phase precipitation fields derived from the Darwin 
volume scan are shown in Figs. 3a and 3b, respectively. Fig. 3a may be compared to the 


SSM/I retrieval of integrated liquid precipitation in Fig. 2a. Precipitation amounts in the 
convective leading edge and trailing stratiform regions of the tropical squall line have 
roughly the same magnitudes in both the SSM/I and radar estimates. A small positive bias 
in liquid precipitation estimates is evident in the SSM/I retrievals, however. Also, there 
appear to be spurious retrievals of light preciptation north of the squall line and along the 
lower edge of Melville Island. The high bias and spurious retrievals may be due, in part, to 
the lack of explicit nonprecipitating clouds in the cloud ensemble/radiative parameterization. 
The effect of nonprecipitating clouds is imbedded in the relationships between precipitation 
amount and upwelling brightness temperature, but the cloud liquid water content does not 
vary as an independent variable. A region containing significant nonprecipitating cloud can 
therefore be interpreted as a region of light precipitation in the retrieval. Spurious 
retrievals of liquid precipitation in coastal regions could also be due to slight errors in 
SSM/I data geolocation, coupled with the very low sensitivity of the SSM/I channels to 
variations in precipitating water in coastal regions. Although the SSM/I geolocation was 
checked and adjusted to match up with a coastline data base, any small error which would 
move the SSM/I data from land over ocean would require a large adjustment in the land 
fraction, fiand» to compensate. An increase in the amount of liquid precipitation in the 
retrieval solution could also provide this compensation. In addition, the derivative scaling 
for liquid precipitation in the vicinity of coastlines must be adjusted to compensate for the 
lack of sensitivity of upwelling brightness temperatures to variations in precipitating liquid 
amount. This is because the derivative goes through a change of sign and actually vanishes 
at a coastal boundary, and therefore derivative magnitudes are relatively small. The scaling 
of derivatives heightens the sensitivity of the parameterized brightness temperatures to 
changes in precipitating liquid, and leads to exaggeration of errors in the liquid precipitation 
retrievals arising from noise or improper geolocation. 


The radar-derived slant path-integrated precipitating ice in Fig. 3b may be compared 
to the SSM/I retrieval in Fig. 2b. There is a fairly good correspondence between the radar 
and SSM/I retrieved Fields, both in pattern and magnitude of precipitating ice. The figures 
suggest only a slight high bias of the SSM/I retrieved ice amounts over the radar estimates. 
Note that within about 20 km of the radar site there is a conical data-void region in the 
radar-derived field, due to the 25° limiting elevation angle of the radar (Keenan, et al., 
1988). 

3.3 Tropical Storm Emily - 1004 UTC Overpass 

The physical retrieval method is applied to SSM/I observations of the incipient 
Tropical Storm Emily at 1004 UTC on September 21, 1987, as it was developing over the 
Windward Islands off the northeast coast of South America. The SSM/I 37 GHz and 85.5 
GHz horizontal polarization images of Emily on this date are shown in Fig. 4a and 4b, 
respectively. The northern coast of South American and Caribbean Islands are indicated by 
white lines, and a brightness temperature grayscale is displayed at the bottom of each panel. 
The panel dimensions are approximately 1400 km x 1400 km. The most intense convective 
development in Emily is seen just north of the island of Trinidad and Tobago. This oblong 
convective zone spans about 250 km, with 37 GHz horizontal polarization antenna 
temperatures near 255 K, indicating strong absorption and emission from precipitating 
liquid. Embedded in the convective zone are smaller regions of microwave scattering at 
85.5 GHz, especially to the north-northeast of Trinidad and Tobago. The more intense 
scattering regions are associated with precipitating ice held aloft by the stronger updrafts of 
the convective zone, and the 85.5 GHz horizontal polarization antenna temperatures reach a 
minimum of 150 K. To the north and northwest of the convective zone are semi-circular 
bands in the antenna temperature imagery arising from developing bands of precipitation. 
Most of these bands contain significant amounts of precipitating ice, identified by scattering 
depressions in the 85.5 GHz antenna temperatures. 


The retrieval algorithm is applied to a 64 x 64 swath grid section centered on Emily. 
After 21 iterations the error functional in Eq. (1) is minimized to yield the retrievals of grid 
box average, slant-path integrated precipitating liquid and ice in Figs. 5a and 5b. It may be 
noted that the pattern of precipitating liquid in Fig. 5a closely follows the pattern of 37 GHz 
emission in Fig. 4a. Maximum values of integrated liquid occur in the convective core area 
of Emily, with magnitudes close to 2.6 kg/m^, and lesser amounts, on the order of 1 
kg/m^, in developing rainbands. Retrievals of precipitating liquid along the coast of South 
America are spurious, for the most part, since coincident infrared imagery (not shown) at 
the time of the SSM/I overpass does not show significant cloud along the coast. The 
precipitation artifacts are likely due to a compensation in the retrieval for variations in land 
surface emissivity. 

Integrated precipitating ice concentrations in the retrieval of Fig. 5b coincide with 
regions of strong scattering in the 85.5 GHz SSM/I imagery. Precipitating ice retrievals 
exhibit maximum values close to 3 kg/m^ in the convective core region, with values up to 1 
kg/m^ in the outer bands. The overall extent and magnitude of the precipitating ice 
amounts are generally less than those of precipitating liquid in Emily at the overpass time. 

3.4 Tropical Storm Emily- 2242 UTC overpass 

The SSM/I retrieval procedure is repeated for a subsequent overpass of Emily at 
approximately 2242 UTC on September 21, 1987. SSM/I imagery from this overpass in 
the 37 GHz and 85.5 GHz horizontal polarization channels is shown in Figs. 6a and 6b, 
respectively. These may be compared to the imagery of Emily roughly 12 hours earlier, 
shown in Figs. 4a and 4b. It may be noted from the figures that the core region of Emily 
had moved about 400 km to the west-northwest of its original position in the Windward 
Islands. Also, the storm had taken on the classical structure of a hurricane, with spiral 
bands and an apparent eye which can be seen in the 85.5 GHz imagery. At the time of the 
second SSM/I overpass, Emily was still classified as a tropical storm, but was deepening 


rapidly, and the storm attained hurricane status early the following day (ref. Case and 
Gerrish, 1988). 

Retrievals of slant-path integrated precipitating liquid and ice from the second 
overpass of Emily are diplayed in Figs. 7a and 7b, respectively. A correspondence 
between precipitating liquid and 37 GHz emission patterns (Figs. 6a and 7a) is evident. 
Also there is a strong correlation between scattering depressions in brightness temperature 
and retrieved precipitating ice amounts (Figs. 6b and 7b), as it was noted in the 1004 UTC 
Emily retrievals. Maximum retrieved values of precipitating liquid and ice within Emily are 
on the order of 2 kg/m^, although liquid precipitation amounts tend to be higher in general. 

4. EXTENSION OF CURRENT WORK 

AND APPLICATIONS 

Work will continue to further improve the representation of clouds and their 
environments in the sensor response model. Greater variability in the solution fields of 
emissivity, air temperature, and water vapor may also be allowed if additional sensor data 
and/or constraints are utilized in the retrieval. Future tests will include the lower-frequency 
SSM/I data as well as SSM/T, SSM/T-2, and geostationary infrared measurements. 
Retrievals will be validated using volume-scan radar measurements and conventional 
weather data. 

The sequence of Emily retrievals is utilized in the latent heating analysis and 
numerical weather prediction model forecasts in Section 4 of this report. Rates of latent 
heating are inferred from the liquid and ice precipitation distributions and are assimilated 
into the forecast model fields using a new technique. This application is important because 
the development or "spin-up" of storms in numerical model forecasts is often retarded due 
to a lack of fine-scale initial data. The SSM/I-derived latent heating information can 
partially compensate for this lack of fine- scale data. 

Latent heating information may help to improve short-term numerical weather 
prediction model forecasts of storms, but forecast model fields may also be used to help 


improve retrievals through the constraint term in Eq. (1). The interplay between forecast 
models and retrievals is a subject of further study by the author. 
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‘ , ^ e . tr ^ eve< ^ geophysical parameters (P) and associated spatial constraints in the 
SSm/I physical retrieval method. 


Parameter 

Spatial Constraints 

land fraction within each 
grid box 

none 

land emissivity parameters 
from Grody (1988) 
parameterization 

uniform over grid 

ocean emissivity parameters 
from Grody (1988) 
parameterization 

uniform over grid 

first eigenfunction weight of 
water vapor mixing ratio 
profile 

uniform over land, 
no constraint over ocean 

precipitation fraction within 
each grid box 

none 

mean slant path depth of 
precipitating liquid within the 
precipitation fraction 

none 

mean slant path depth of 
precipitating ice within the 
precipitation fraction 

none 



(a) 




Fig. 1. SSM/I antenna temperature imagery of a tropical squall line in the vicinity of 
Darwin, Australia at 0906 UTC on February 9, 1988. The 37 GHz horizontal polarization 
imagery is shown in panel (a), and the 85.5 GHz horizontal polarization imagery is shown 
in panel (b). Coastlines are indicated by solid white lines. The area represented by each 
panel is approximately 250 km x 250 km. 





Fig. 2. Retrieved distributions of slant path-integrated precipitating liquid (a) and ice (b) 
derived from the SSM/I overpass of Darwin, Australia at 0906 UTC on February 9, 1988. 




Fig. 3. Distributions of slant path-integrated precipitating liquid (a) and ice (b) derived 
from a volume scan of the Darwin/TOGA radar starting at 0850 UTC on February 9, 1988. 




Fig. 4. SSM/I antenna temperature imagery of Tropical Storm Emily in the Windward 
Islands at 1004 UTC on September 21, 1987. The 37 GHz horizontal polarization imagery 
is shown in panel (a), and the 85.5 GHz horizontal polarization imagery is shown in panel 
(b). Coastlines are indicated by solid white lines. The area represented by each panel is 
approximately 1500 km x 1500 km. 




(a) 




Fig. 5. Retrieved distributions of slant path-integrated precipitating liquid (a) and ice (b) 
derived from the SSM/I overpass of Tropical Storm Emily at 1004 UTC on September 21, 
1987. 



5g- 6 - SSM/I antenna temperature imagery of Tropical Storm Emily north of Venezuela at 
2242 UTC on September 21, 1987. The 37 GHz horizontal polarization imagery is shown 
m panel (a), and the 85.5 GHz horizontal polarization imagery is shown in panel lb). 
Coastlines are indicated by solid white lines. The area represented by each panel is 
approximately 1500 km x 1500 km. p 




Fig. 7. Retrieved distributions of slant path-integrated precipitating liquid (a) and ice (b) 
derived from the SSM/I overpass of Tropical Storm Emily at 2242 UTC on September 21, 
1987. 
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Abstract 


Retrieved three dimensional liquid water fields are used in a regional mesoscale forecast model 
to help induce storm development and enhance the early production of precipitation. Vertical layer- 
mean heating rates are regressed against slant path-integrated liquid and ice precipitating water 
contents, based upon a cloud ensemble model. This procedure yields a best fit two-parameter 
regression formula for each vertical layer. The regression formulae are used to estimate latent 
heating values from the slant-path integrated liquid and ice precipitation water contents from SSM/I 
retrievals. In the forecast model, diabatic temperature contributions are calculated and used in a 
diabatic forcing procedure with or without diabatic initialization. We find that the time needed to 
spin-up precipitation processes in tropical storm Emily is greatly accelerated, and the horizontal 
distribution of the forecast fields is improved through the application of the data. 
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1. Introduction 

Mesoscale forecasts must go through a development stage when initial conditions do not 
adequately resolve mesoscale circulations. This spin-up is common in numerical weather prediction 
since much of the data is obtained from synoptic scale observations. Also, the failure to include 
cloud, rain and ice fields in the assimilation procedure means than this spin-up process is repeated 
with the next sequential forecast This neglect is a natural consequence since cloud and rain water 
data have been unavailable. There has been some effort to circumvent this problem through the use 
of diabatic or cumulus initialization techniques (Errico and Rasch 1988; Donner 1988; Turpeinen et 
al. 1990; Tarbell et al 1981; Salmon and Warner 1986, etc.). These procedures have been 
somewhat successful, but the vertical distribution of the diabatic heating remains a major problem 
since raingauge reports give a measure of surface precipitation depth, and visible/infrared remote 
sensing methods typically yield estimates of surface rain rate. 

Some success in assimilating explicit cloud and rain water information, and initializing diabatic 
processes was demonstrated in chapter one of this report. A continuation of these procedures is 
presented here, but using radiance data measured by the Special Sensor Microwave/Imager (SSM/I) 
to estimate the heating rates. SSM/I-derived estimates of integrated precipitating liquid and ice are 
optimized using physical information supplied by a cloud ensemble/radiative model of convective 
and stratiform clouds. Vertical layer-mean latent heating rates from the cloud ensemble model are 
regressed against the slant path-integrated liquid and ice precipitation water contents to determine the 
best fit two parameter regression coefficients for each vertical layer. Given SSM/I retrievals of path 
integrated liquid and ice precipitation amounts, the regression formulae yields a vertical distribution 
of heating rates for our forecast model applications. 

In this final report only the data associated with tropical storm Emily is evaluated. This will be 
expanded in the near future for a proposed journal article. Forecasts are made using the National 
Meteorological Centre’s (NMC) initialized global analyses. When only an adiabatic initialization is 
used, our forecasts produce no development in the vicinity of the Emily disturbance. This result 
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suggests that the initial data are missing most of the mesoscale information needed to begin the storm 
development Forecasts are dramatically different when either diabatic forcing or diabatic forcing 
combined with diabatic initialization is used. These findings are presented and discussed in this 
portion of the final report but first, additional details about the three dimensional diabatic heating rate 
data are presented. 

2. Latent Heating Rates for Emily at 1004 UTC, September 21, 1987 

Latent heating rates in Emily at 1004 UTC, September 21, 1987 are derived from the SSM/I 
retrieval of the slant path-integrated precipitation liquid and ice distributions discussed in chapter 3 of 
this report. The latent heating analysis is based upon regression formulae which are applied to the 
retrieved precipitation fields. The regression formulae are derived as follows. 

Latent heating rates from the convective and stratiform cloud simulations in each of the cloud 
ensembles of chapter 2 are interpolated to grid points within the 25 km x 25 km x 20 km domain 
shown in Fig 3 of that section. The heating rates are then averaged in 1.1 km thick layers over the 
25 km x 25 km domain to obtain layer-mean estimates at 17 levels for each cloud ensemble. The 
layer-mean latent heating rates are regressed against the slant path-integrated liquid and ice 
precipitation water contents in the same ensemble using the formula 

Hiat =a 1 *PLWP -l- a 2 *PIWP , (1) 

where H lal is the latent heating rate in joules/m 3 -s, PLWP is the slant path-integrated liquid 
precipitation averaged over the grid box, PIWP is the slant path-integrated ice precipitation averaged 
over the grid box, and aj and a 2 are the best-fit regression coefficients. Both the slant path- 
integrated liquid and ice amounts are assumed to have units of kgnr 2 in (1). PLWP and PIWP are 
chosen as predictors of the latent heating rate, since these parameters are estimated in the SSM/I 
physical retrieval method (see chapter 3). Note that H lat in the above formula equals zero if both 


4 


PLWP and PIWP are zero, which is a constraint on the regression formula since it is assumed that 
their is no latent heating if no cloud or precipitation exist. 

A total of 213 convective cloud ensembles and 42 stratiform cloud ensembles based upon the 
GATE day 261 environment (1200 UTC sounding from Quadra research vessel) are utilized in the 
regression analysis. The regression coefficients and formula errors for each of the 17 standard 
levels are summarized in Table 1. Note that at levels below the freezing level (4km), latent heating 
rates increase with the path-integrated precipitation liquid amount and decrease with the path- 
integrated precipitating ice amount The negative correlation between latent heating rates and ice 
amount (negative values of a 2 ) can be linked to the differences between the thermodynamics of 
convection and stratiform precipitation. Low level heating is normally associated with convective 
precipitation events (Tao et al. 1990). However, if precipitating ice concentrations are high relative 
to precipitating liquid amounts, then a stratiform precipitation event is indicated. Stratiform 
precipitation is usually associated with cooling at altitudes below the freezing level, due to melting of 
precipitating ice and evaporation of rain. 

At levels between 5.70 km and 10.70 km, the latent heating rate is positively correlated with 
both the liquid and ice precipitation amounts. At higher altitudes the magnitude of the latent heating 
rate fall off substantially, and it is difficult to infer any systematic trends. Latent heating rates at 
17.45 km and 18.58 km are so low that meaningful regressions are not possible. It may also be 
noted that correlations between latent heating rate regression estimates and the true heating rates are 
significant at and below 10.70 km, with the exception of the 2.82 km and 3.94 km levels. 

The latent heating rates are inferred from the slant-path integrated precipitating liquid and ice 
estimates obtained from the SSM/I physical retrieval method, upon inserting these estimates into the 
regression formula (1) above. The latent heating rates at the standard levels (see Table 1) are then 
interpolated to levels at 50 mb increments between 1000 mb and 100 mb inclusive. These 
vertically-interpolated values are then interpolated in the horizontal from the 25 km resolution swath 
grid to the 40 km resolution forecast model grid. The contour plot in Fig 1 , superimposed on the 
GOES IR image valid at 1001 UTC 21 September 1987, shows the distribution of temperature 


5 


increase associated with latent heating at sigma level 10, approximately the 700 mb level, inferred 
from the 1004 UTC Emily retrieval (also see Figs 5a and 5b of chapter 3). Here the maximum latent 
heating produces a temperature change of 0.07 C per 300 s' 1 . The bi-modal vertical distribution of 
the latent heating is largest at sigma levels 12 and 6, approximately at the 850 and 500 mb levels, 
where the maximum diabatic activity is equivalent to 0.44 C and 0.37 C per 300 s' 1 , respectivity. 

3. The Diabatic Procedure 

A diabatic initialization and/or diabatic forcing procedure is installed into the CIMSS forecast 
model to take full advantage of the semi-implicit method, as previously described in chapter 1 . This 
is accomplished by modifying the semi-implicit average temperature T 2t , computed from the 
thermodynamic equation (McGregor et al 1978; Leslie et al 1985), by adding the externally 
determined temperature contribution T ex . This redefines the average temperature at the appropriate 
grid points according to 


T 2t _ X 2t + Xm f (2) 

where 

T ex =(LAt/pc p )H lat , (3) 

and 

X 2x = [X t + i + -p-ij/2 . (4) 

Here c p is the specific heat for constant pressure Ji lat the remotely sensed heating rate determined by 
(1), L is the latent heat of condensation or evaporation and the superscripts signify the time step. 
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During the diabatic forcing or diabatic initialization process, which is designed to initialize the 
mesoscale circulation, the mixing ratio can be altered as desired To avoid excessive evaporation, 
some upward adjustment of the mixing ratio may be required, namely 

q v = q v + Aq v . (5) 

In our calculations Aq v is set to give the equivalent of the amount of moisture needed in the 
condensation process to achieve the observed heating rates. 

During the vertical mode initialization procedure (Bourke and McGregor 1983) T ex is held 
constant at the selected grid points. The application early in the forecast period of an external 
diabatic forcing represents a natural extension of the diabatic initialization. This procedure requires 
at each time step that the grid point average temperature is modified by the difference between the 
external contribution T ex and the model produced diabatic temperature T IN = |AtLQ/c p I, yielding 

T 2t a T 2t + (T„ - Tjn) . (6) 

Here Q is the model produced liquid water condensation at each grid point over one time step. In 
our study the quantity in parenthesis in (6) is taken to be non-negative. Equation (6) can be 
interpreted as a blending over time, between the forcing and the model diabatic production. If all 
precipitating grid locations have sufficient condensation, then the external forcing is reduced to zero. 
Even though this is the desired result it is not explicitly required before the termination of the forcing 
at a preselected forecast time. This application of the external forcing is somewhat modified from the 
procedure used in Newtonian nudging (Hoke and Anthes 1976). 

Note that (2) and (6) are similar, except model produced condensation reduces the external 
forcing in (6). This reduction removes the potential for over forcing. The external diabatic forcing 
is a natural continuation of diabatic contributions introduced in the initialization. The forcing is 
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applied starting with the first time step after the vertical mode initialization and continued until the 
model produces an equivalent temperature change or until a preselected termination time is exceeded. 
With each passing time step, T ex is advected by the model horizontal winds according to 

a(T„)/dt + u3(T ex )/3x + v0(T ex )/5y = 0. (7) 

In regions with strong vertical shear it may be advantageous to advect T ex using a vertically averaged 
horizontal wind. The advection process will ultimately smooth and diffuse the heating if (7) is used 
over a sufficiently large time period. Normally the diabatic forcing lasts only a few hours. In any 
event, some tapering must be used to terminate the diabatic forcing, provided (6) is not reduced to 
zero by model generated condensation processes. Remotely sensed time variations in T ex could 
also be allocated for if known. In our calculations, T ex is assumed to remain unchanged except for 
horizontal advection. 

An alternative approach is to apply the external diabatic forcing at the satellite observation time. 
Then, a gradual transition into and exit from the external diabatic forcing regime gready avoids the 
production of noise and gravity waves during the forecast. The time span needed to produce 
gradual changes depends upon the magnitude of the external forcing. The length of time period 
used in the forcing can be varied to achieve the desired result. Excessively long periods of forcing 
should be avoided since they can can lead to some negative consequences. For example, too much 
diabatic forcing can cause a build up of excessive amounts of liquid water that can eventually lead to 
large amounts of precipitation drag which cause downbursts that choke off the convergence and the 
moisture supply. Events of this type can reduce the overall storm development because they may 
produce storm-wide low level divergence that chokes off the moisture supply. 

4. Forecasting Emily 

In Figs 2 and 3 are 700 mb geopotential heights and divergence fields from a twenty four hour 
forecast beginning at 0000 UTC 21 September 1987. These fields are super-imposed over the 
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GOES infrared (IR) image for the region containing tropical storm Emily valid at 0001 UTC 22 
September 1987. This simulation uses an adiabatic initialization with 4 vertical modes. Note in 
Figs 2 and 3 that the height and divergence fields show no development whatsoever in the immediate 
vicinity of Emily. No precipitation is associated with this flow. However, some precipitation does 
occur over elevated topographical features in association with changes in the diurnal cycle. 

Retrieved heating rates for Emily, based on the physical retrieval and regression procedures 
described in section 2, provides a knowledge of both the horizontal and vertical distribution of the 
diabatic heating. This information is then used directly in the CIMSS model diabatic initialization 
and diabatic forcing procedures, as described in section 3. Diabatic forcing, beginning at time step 
120, i.e., at 1000 UTC 21 September 1987, and ending at time step 165, is applied during a forecast 
according to (2) through (7). The introduction contains a one and one-half hour ramp-up, with one 
and one-quarter hours of constant forcing, and a one hour ramp-down procedure. Thus, the diabatic 
forcing begins slowly ten hours into the forecast and shuts down 3.75 hours later. The length of the 
ramps can be adjusted based on model response. A tapered introduction and exit is necessary to 
keep gravity and inertia wave noise to a minimum. 

At the end of the twenty-four hour forecast, or fourteen hours after the start of the forcing, the 
700 mb geopotential heights contain several closed contour centers, as seen in Fig 4. Some 
indication of the development of a warm core is provided by the center containing the 3200 m 
contours. The convergence (negative divergence) at 700 mb (Fig 5) has a maximum of 
approximately -8 x 10' 5 s _1 directly over the location of Emily. Also, the positive vorticity reaches 
a magnitude of 24 x 10 5 s 1 at 700 mb, as illustrated in Fig 6. The integrated precipitating liquid 
water, fourteen hours after the start of the diabatic forcing, has a maximum of 4.0 kgm* 2 (Fig 7). 
However, smaller magnitudes are found in the SSM/I physical retrieval, shown in Fig 8, which 
contains a maximum near 2.5 kgm 2 of precipitating liquid and ice. Also, the spatial distribution for 
the retrieved precipitating field shows a comma shape. The lack of spatial structure in the forecasted 
quantity is though to be linked to a common feature among models, namely, a tendency to over- 
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develop precipitation situations. Too much condensation occurs in grid scale explicit cloud 
calculations because moist and dry downdrafts, which help to restrict the total condensation, are not 
correctly parameterized. Inadequacies in the vertical and horizontal mixing or turbulence 
parameterization, too coarse of a grid resolution, and other deficiencies in the cloud and rain water 
parameterization s may also be contributing factors. To help keep the storm from over developing or 
developing too rapidly, in our calculations the maximum amount condensed at any one time step at 
any grid location is restricted by requiring that 0< lQ/fq l<0.5, where f is the Coriolis parameter and 
q the mixing ratio. Conservation of moisture is retained since the excess is carried in the slightly 
super saturated mixing ratio. This is thought to be a better alternative than using artificial diffusion, 
which is non-conservative, in our 40 km resolution simulations. 

The remotely retrieved heating rate information described in section 2 is used in a fourteen hour 
forecast containing a diabatic initialization followed by two hours (24 time steps) of diabatic forcing. 
Four vertical modes are used in this simulation since, as described in our earlier work, the combined 
diabatic initialization and diabatic forcing is not sensitive to the choice of the vertical wave number in 
the vertical mode initialization (Bourke and McGregor 1983). 

Results from the fourteen hour forecast are shown in Figs 9 through 12. The forecasted 700 
mb height field contains two centers with closed contours (Fig 9). For the same time period, 
convergence (-1.2 x 10' 5 s 1 ) is found over this region (Fig 10) associated with a positive vorticity 
that has a maximum value of 16 x 1 O' 5 s 1 (Fig 11). This development, that coincides with Emily, is 
somewhat less than that shown in Figs 4 through 8 primarily because the length of the forcing is 
less. This weaker signal is also evident in the integrated rain water shown in Fig 12, when 
compared with the values shown in Fig 7. Dispite having less integrated rain water, the weaker 
system has more rainfall at the surface during the first few hours of its existence, presumely because 
the clouds are closer to the surface and because the weaker vertical motions do not gready impede the 
precipitating liquid. 
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5 . Summary 

Vertical heating profiles (1000 UTC 21 September 1987) are constructed using the vertical 
layer-mean latent heating rates from a cloud ensemble model and the slant path-integrated liquid and 
ice precipitating water contents retrieved from the SSM/I. This information is then used in a diabatic 
initialization and/or diabatic forcing procedure to help encourage the spin-up the precipitation and the 
development of tropical storm Emily. The robust nature of the diabatic forcing procedure, as 
described and tested in chapter 1 of this final report, is evident by the successful generation of a 
disturbance in both the twenty four and fourteen hour forecasts. One additional case study should 
be done to complete our investigation, after which the diabatic forcing procedure should be tested in 
an operational setting where many different meteorological situations can occur. The physical 
retrieval technique for estimating the vertical structure of the latent heating looks extremely 
promising. 
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Table 1. Latent heating rate regression formulae coefficients and statistics. Note that there are 
insufficient data to perform a regression at the two highest levels. Units of altitude are 
kilometers, and bias and are in joules/m 3 -s. 


altitude 

at 

a 2 

bias 


r 

0.56 

0.370 

-0.251 

-0.073 

0.430 

0.617 

1.69 

0.909 

-0.355 

-0.087 

0.942 

0.528 

2.82 

0.108 

-0.122 

0.066 

0.478 

0.229 

3.94 

0.227 

-0.119 

0.018 

0.532 

0.305 

5.07 

0.344 

0.044 

0.007 

0.289 

0.607 

6.19 

0.166 

0.168 

0.008 

0.185 

0.786 

7.32 

0.045 

0.181 

0.007 

0.139 

0.841 

8.45 

0.059 

0.123 

0.007 

0.112 

0.809 

9.57 

0.044 

0.072 

0.009 

0.062 

0.853 

10.70 

0.038 

0.025 

0.006 

0.054 

0.650 

11.82 

0.045 

-0.003 

0.006 

0.065 

0.436 

12.95 

0.170 

-0.022 

0.022 

0.162 

0.576 

14.08 

0.041 

-0.007 

0.002 

0.049 

0.461 

15.20 

-0.002 

0.000 

- 0.000 

0.002 

0.545 

16.33 

17.45 

18.58 

- 0.000 

0.000 

0.000 

0.000 

0.443 
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List of Figures 


Fig. 1 . Temperature values (C) x 100 used in the diabatic initialization and/or diabatic forcing 
routine at sigma level 8 (approximately 700 mb). 

Fig. 2. Height field at 700 mb from the twenty four hour forecast with adiabatic initialization 
conditions and no diabatic forcing. Contour interval is 4 m. 

Fig. 3. Same as Fig 2 but for the divergence field. Contour interval is 1 x 10' 5 s 1 . 

Fig. 4. The 700 mb height field from the twenty four hour forecast with diabatic forcing applied 
beginning 10 hours into the forecast and lasting for 3.75 hours. Contour interval is 4 
m. 


Fig. 5. Same as Fig 4 but showing the divergence field. The contour interval is 1 x 10 5 s 1 . 

Fig. 6. Same as Fig 4 but showing the vorticity. Contour interval is 4 x 10' 5 s 1 . 

Fig. 7. Vertically integrated rain water x 10 associated with the fields shown in Figs 4, 5 and 6. 
Contour interval is 0.5 kgm" 1 . 

Fig. 8. Retrieved estimates of the integrated precipitating liquid and ice from SSM/I radiances 
(Olson, 1993). 

Fig. 9 The 700 mb heights from the fourteen hour forecast that used diabatic initialization and 
two hours of diabatic forcing. Contour interval is 4 m. 
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Fig. 10. Same as Fig 9 but showing the divergence. Contour interval is 1 x 10' 5 s 1 . 


Fig. 1 1. Same as Fig 9 but showing the vorticity. Contour interval is 4 x 10' 5 s 1 . 


Fig. 12. The vertically integrated rain water x 10 associated with the fields shown in Fig 9, 10 and 
11. Contour interval is 0.5 kgm 1 . 


16 




























Reprinted from Monthly Weather Review, Vol. 121, No. 7, July 1993 

American Meteorological Society 


N 9 


Moist Wind Relationships 


0 7 

/ i /3#3 


William H. Raymond 

Cooperative Institute for Meteorological Satellite Studies, University of Wisconsin, Madison, Wisconsin 
(Manuscript received 24 July 1992, in final form 8 January 1993) 


ABSTRACT 

Equations describing the temporal and spatial behavior of the kinematic moisture and heat flux are introduced 
in this study. In these nonlinear equations, the contribution by diabatic processes to the large-scale flux is 
composed of two parts. One part is associated with a Rayleigh damping term, while the other arises from 
temporal and spatial changes in the pressure gradient term. 

The influence of diabatic processes on the large-scale moisture fluxes depends greatly on the degree of balance 
between forcing and damping terms in the governing equations. The existence of a near balance requires a 
reduction in the large-scale horizontal geostrophic wind speed. From a scale analysis of the moisture flux equations 
it is argued that reductions in the large-scale horizontal wind speed, observed within major cumulus cloud 
systems, help conserve large-scale moisture fluxes. The deviation of the wind from geostrophic conditions is 
easily estimated. This wind modification induces secondary vertical circulations that contribute to the convergence, 
creating or supporting long-lived mesoscale flows. In the tropics the wind modification has an antitriptic rela- 
tionship. 

These diagnostic findings suggest possible modifications to the wind field in the application of cumulus 
parameterization, and may be important in diabatic initialization of numerical weather prediction models. 


1. Introduction 


In the extratropics, the large-scale atmospheric flow 
fluctuates about known balance conditions. Away from 
the major centers of cyclonic and anticyclonic activity, 
the horizontal wind components closely approximate 
geostrophy in the free atmosphere. Within strong cy- 
clonic centers the flow is modified by curvature and 
by diabatic effects. To forecast the wind field in the 
presence of diabatic activity would require solutions 
of the complete equations of motion with physical pro- 
cesses included, which is not always practical. Other- 
wise, to get a snapshot of the instantaneous horizontal 
wind field requires solving the nonlinear diabatic bal- 
ance equations. The simple wind laws, like the gradient- 
wind approximation, do include the effect of curvature 
and some indirect consequences of diabatic heating. 
This study outlines a simplified approach derived from 
first principles that identifies how diabatically induced 
modifications in the horizontal wind components in- 
fluence the moisture and sensible heat fluxes. This 
study provides necessary additional information about 
the effect of diabatic processes on the atmospheric dy- 
namics for use in model initialization and cumulus 
parameterization . 

Many years of observations have shown that diabatic 
processes can significantly enhance cyclonic develop- 
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ment and influence the general circulation. Many 
journal articles and textbooks testify to this. A small 
subset includes Petterssen ( 1956), Riehl and Malkus 
(1958), Kuo (1965, 1974), Yanai et al. ( 1973), Riehl 
(1979), Donner et ai. (1982), Houze and Hobbs 
(1982), and Konig and Ruprecht (1989). Modeling 
studies also have verified that moist processes are im- 
portant for rapid development (Kuo and Reed 1988; 
Reed et al. 1988). Additionally, diabatic initialization 
techniques are used to enhance cyclonic circulations 
and low-level convergence to help speed up precipi- 
tation processes in numerical forecast models (e.g,, 
Tarbell et al. 1981; Krishnamurti et al. 1983; Krish- 
namurti et al. 1984; Donner 1988; Wang and Warner 
1988; Turpeinen et al. 1990; Daley 1991 ). 

There also are many studies on the effect of clouds 
on the large-scale dynamics, including Malkus ( 1952), 
Schneider and Lindzen ( 1976), Sui and Yanai ( 1986), 
Cho (1985), and Konig and Ruprecht (1989). The 
vertical transport of atmospheric properties in clouds 
is important because the cloud vertical velocity is much 
larger than that found in the environment, especially 
in the tropics. In contrast, the magnitudes of the hor- 
izontal flow in the cloud and environment are of the 
same order, but according to the above studies, mixing, 
vertical transport, and other processes are responsible 
for some deceleration of the mean flow by “cumulus 
friction” (Malkus 1952; Schneider and Lindzen 1976). 
Several articles propose parameterizations of cumulus 
convection effects on the large-scale vorticity or mo- 
mentum. Approaches include mass entrainment and 
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detrainment concepts (Ooyama 1971 ), vertical trans- 
port and properties from simple cloud models (Shapiro 
and Stevens 1980), the introduction of drag terms (Es- 
bensen et al. 1987; Austin and Houze 1973), and other 
representations of cumulus friction (Schneider and 
Lindzen 1 976 ) . In this study a different method is pro- 
posed to estimate the cumulus friction process. 

Gill (1980) has examined tropical circulation forced 
by localized heating, while Heckley and Gill (1984) 
have studied transient features in the tropics when the 
heating is impulsively started. In these studies the role 
of diabatic heating on the large-scale flow is modeled 
by adding a forcing term to the continuity equation in 
the linearized shallow-water model. The forcing mod- 
ifies the pressure and its gradients, which in turn modify 
the horizontal wind field. This model is conceptually 
simple yet captures the essence of the large-scale re- 
sponse to diabatic forcing arising from thermal changes. 
Gill (1980) heuristically added a Rayleigh damping 
term with constant coefficient over the entire domain 
of his linear model, which produced a steady state, but 
this was not intended to represent the cumulus friction 
process, and it provided no insight on why cumulus 
friction occurs, its region of influence, or consequences. 
Holton and Colton (1972) in their linear vorticity 
model also found that damping was needed to obtain 
agreement with observations. 

Lorenz (1978) has identified some optimal atmo- 
spheric energetic arrangements, but local optimal re- 
lationships remain unknown if they exist at all. One 
can speculate that cumulus friction is partially an at- 
mospheric response for conserving basic large-scale at- 
mospheric properties or optimizing fundamental pa- 
rameters like the kinetic energy, enthalpy, etc. Clearly, 
temporal and advective changes in moisture flux, heat 
flux, and entropy are influenced by cumulus friction. 

In this diagnostic study it is found that mesoscale 
moist processes introduce into the moisture and heat 
flux equations a Rayleigh damping term that depends 
upon the condensation and evaporation rate. This fric- 
tion allows changes in the large-scale moisture fluxes 
to remain relatively small provided a near balance exists 
between the forcing and damping terms. Scale analysis 
shows that this balance involves the pressure gradient, 
Coriolis, and Rayleigh damping term. Balance further 
requires a reduction in the large-scale horizontal geo- 
strophic wind speed, that is, the creation of an ageo- 
strophic wind. Clearly, if the horizontal wind is reduced 
excessively, then the balance condition breaks down. 
Our hypothesis is that for long-lived mesoscale con- 
vection the large-scale horizontal winds do not differ 
greatly from that needed to establish a near balance 
between the pressure gradient, Coriolis, and Rayleigh 
damping terms. It is found diagnostically and in nu- 
merical applications that the ageostrophic wind en- 
hances the existing convergence and establishes a me- 
soscale vertical circulation within the large-scale flow 
to help balance the modified horizontal wind. 


To illustrate our hypothesis we examine the behavior 
of the kinematic moisture and heat flux. The derivation 
of these equations is quite simple, but they yield sur- 
prising information. In these equations the contribu- 
tion of diabatic heating to the horizontal wind com- 
ponents is naturally partitioned into two parts. One 
diabatic contribution is contained within the geo- 
strophic wind component (pressure gradient term), 
while the remaining one is associated with the Rayleigh 
damping (condensation ) or forcing (evaporation ) term. 
To examine the importance of the terms in the flux 
component equations, a scale analysis is preformed and 
simplified solutions are found. For nontropical regions 
our approximate solutions are given in terms of the 
geostrophic flow. In some of our illustrative examples 
a time scale less than the inverse of the Coriolis param- 
eter is required; however, the flux equations are valid 
even when this is not true. In tropical regions any ten- 
dency toward moisture flux conservation enhances the 
relationship between the pressure gradient term and 
the Rayleigh damping term. This induces secondary 
antitriptic circulations. A simple hurricane model is 
presented to illustrate the wind modifications associated 
with any enhancement toward moisture flux conser- 
vation. Also tested are applications in cumulus param- 
eterization. 


2. The equations 


The inviscid equations (Haltiner and Williams 1980, 
pp. 16, 17, and 308 ) for large-scale flow describing hor- 
izontal wind components u and v y the temperature T, 
potential temperature 0, the mixing ratio q , and pres- 
sure p are 


du 

~Jt 

dv 

dt 

dd 

dt 


1 dp _ 

T" + fo* 

p dx 

1 dp f 
P dy 


L8Q 
Tc p ’ 



( 1 ) 

( 2 ) 

(3) 

(4) 


The total derivative d( )/dt , x, y, z, and t have the 
standard definitions. The Coriolis parameter / satisfies 
/ = 2Q sin0, where ft is the angular velocity of the 
earth and <f> is the latitude. Here p is the density, g is 
the acceleration due to gravity, c p is the specific heat 
for constant pressure, L is the latent heat of conden- 
sation or evaporation, and Q is the condensation or 
evaporation rate. Here Q is a sink for condensation 
and a source for evaporation and sublimation from 
saturated air; otherwise, Q = 0 when q is below satu- 
ration. This is expressed formally by 



July 1993 


RAYMOND 


1994 


<2 = 


dq, 

dt ' 


( 5 ) 


- \_dp 

ffidy- 
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Evaporation from the earth’s surface is not explicitly 
included in the foregoing inviscid model. Adding the 
continuity equation to the four preceding equations 
enables the vertical motion to also be computed. For 
large-scale flow the hydrostatic and ideal gas laws com- 
plete the list of equations. 

The transfer of moisture or heat per unit area per 
unit time represents the moisture or heat flux. The 
kinematic flux form is consistent with atmospheric 
measurements of temperature, mixing ratio, and wind 
speed and direction (Stull 1988). For example, the ki- 
nematic moisture flux has components qu, qv, and 
qw, while the same quantities for the kinematic sensible 
heat flux are 8u, 8v , and 8w. To find the component 
equations governing the temporal and spatial distri- 
bution of the horizontal kinematic moisture flux, mul- 
tiply ( 1 ) and (2) by the mixing ratio q and rearrange 
by application of the chain rule and (4). (A similar 
process involving the potential temperature yields the 
kinematic heat flux equations described in the Appen- 
dix.) The equations for the horizontal components of 
the total moisture flux are 


djqu) 

dt 


q dp 

- -f + qfi + Aqu, 

p OX 


djqi i) 
dt 


-y- qfu + Aqv. 
P 


( 6 ) 

( 7 ) 


Along with the total derivative, note that each nonlinear 
flux component equation contains a pressure gradient 
term, a term containing the Coriolis parameter, and a 
term that includes the factor A = Q/q . The A 
coefficient 1 is negative for condensation, making Aqu 
and Aqv Rayleigh damping terms. Thus, regardless of 
the sign of qu and qv , the damping term always reduces 
the magnitude of the flux. In contrast, evaporation in- 
duces a forcing response that enhances the magnitude 
of the flux. In tropical regions the term containing the 
Coriolis parameter / is small and may be ignored 
without introducing large errors. Then, changes to the 
kinematic flux depend on the residual of the pressure 
gradient and damping terms. 

Away from the tropics the pressure gradient can be 
expressed in terms of the large-scale geostrophic wind 




\_dp 

fpdx" 


( 8 ) 


1 Both the moisture and heat flux have a term containing the con- 
densation or evaporation rate. Rayleigh damping reduces the moisture 
flux, while the heat flux is enhanced by diabatic forcing during pre- 
cipitating events. In magnitude, A is significantly larger than the B 
coefficient associated with the forcing in the sensible heat flux equa- 
tions given in the Appendix. 


It should be noted that diabatic contributions to the 
pressure field are implicitly included in the definition 
of the geostrophic wind. In addition, the Rayleigh 
damping term is important, while other indirect and 
feedback diabatic influences on the large-scale flow oc- 
cur, but are not considered in this study. Using the 
geostrophic approximation simplifies (6) and (7), 
yielding 

d(qu) 

d{ = ~qf(v g - v) + Aqu, (10) 

d(qv) 

-~y- = qf(u g - u) + Aqv. (11) 

The moisture flux is conservative in ( 10) and ( 1 1 ) for 
adiabatic and frictionless conditions, and for geo- 
strophic flow. If the wind components were to remain 
geostrophic during a condensation event, then signif- 
icant reductions in the kinematic moisture flux would 
occur because of the cumulative effects of the Rayleigh 
damping. The exact nature of any optimal local re- 
sponse in the moisture flux, heat flux, energy, or other 
basic atmospheric state is unknown. Scale analysis 
(given below), however, suggests that the large-scale 
moisture flux is nearly conserved. Even though the 
surface evaporation process slowly builds up the at- 
mospheric water vapor, condensation can quickly de- 
plete local supplies unless there exists a mesoscale cir- 
culation to help replenish the local reservoir. Any hy- 
potheses that say there is a tendency toward 
conservation of moisture flux must also include some 
mechanism for generating vertical circulations. 

One way to make the left-hand sides of ( 10) and 
(11) small in magnitude is to alter the wind compo- 
nents u and v so that the right-hand sides are closer to 
being balanced but not necessarily identically balanced. 
It is our hypothesis that the reduction in the large-scale 
horizontal wind components observed during major 
long-lived precipitating events help preserve the large- 
scale moisture flux. This reduction in velocity speed 
decreases the change in the synoptic-scale kinematic 
moisture flux that would otherwise occur if the wind 
speed remained unchanged. As will be shown in the 
example and applications in sections 3 and 4, this 
modification also enhances the longevity of the dis- 
turbance and minimizes loss in the moisture flux by 
creating a mesoscale vertical circulation that helps re- 
store or pump additional moisture into the storm. 
Much of the discussion that follows examines the na- 
ture of the changes in the large-scale horizontal wind 
components that are needed to support our hypothesis, 

A scale analysis of ( 1 ) and (2) shows that for large- 
scale conditions the pressure gradient and Coriolis force 
terms are larger in magnitude than the contribution 
from the total derivative (Holton 1979). In a similar 
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Table 1. A scale analysis* of the terms in the 
u moisture wind [Eq. (8)]. 



c Aqu ) 
dt 

qfv 

qdp 

p dx 

Aqu 

Scaling 

qtf_ 

q2QU 

q Ap 

V 2 q 

~ 

Magnitude 

(extratropics) 

L 

q X 10-* 

q X 10- 3 

P L 

q X 10‘ 3 

<q X 10 -J 

(tropics) 

q X 10' 4 

<q X 10~ 4 

q X I0' 3 

<q X 10~ 3 


* U — 10 m s‘\ A — — s" 1 , L 10 6 m, L c (precipitation scale) 
10 5 m, A p ^ 10 3 N m -2 . 


analysis for (6) and (7) the term containing diabatic 
processes is an order of magnitude more important 
than the contribution from the total derivative, pro- 
vided A is scaled with Q proportional to 10 -7 g g" 1 s” 1 , 
a magnitude common to mesoscale precipitation 
events. Table 1 illustrates our findings for large-scale 
circulations containing a long-lived mesoscale diabatic 
process, of length scale L c , for tropica! and nontropical 
situations. In Table 1 the derivative term yields a mag- 
nitude U 2 q/L , while the diabatic term gives U 2 q/L c , 
where Q is scaled by Uq/L c and L c = 10 5 m. In tropical 
regions the Coriolis parameter is reduced; thus, the 
pressure gradient and the damping terms dominate 
provided A is sufficiently large. This scale analysis sug- 
gests that the large-scale fluid flow is nearly balanced 
by the right-hand sides in (6) and (7). This is consistent 
with our stated hypothesis. Clearly, different scaling 
will provide different magnitudes; for example, short 
duration or microscale diabatic processes are of much 
less importance to the large-scale flow. However, the 
Rayleigh damping term remains important provided 
LJL — 0.1. Convective parameterizations assume a 
similar ratio between grid and diabatic scales. However, 
the total derivative in { 10) and (II) becomes important 
as the scale size is reduced. 

In this study ( 10) and ( 1 1 ) or (6) and (7) are eval- 
uated diagnostically. The right-hand sides of these 
equations are assumed to be nearly in balance, but 
some contribution from the total derivative is allowed. 
A simplified analysis is used to trace how the noncon- 
servative total derivative contributes. As an example, 
assume that the deviation from conservation is known 
to be a small fraction of the flux; then for nonzero 
values of <7, u, and v there exist coefficients a u and a v 
such that 


d(gu) 

dt 

<Hqv) 

dt 


= a u qu. 

(12) 

= a v qv. 

(13) 


This is a generalization of the approximation for the 


partial time derivative used in Gill ( 1 980) and Matsuno 
(1966). To get a grasp of the magnitude of the coef- 
ficient assume that the flux qu experiences a complete 
reversal of sign for the synoptic flow within a time pe- 
riod of A t\ then the coefficient a u would have magni- 
tude 2/ At. Generally the magnitude is less than this. 
The coefficients a u and a v , however, are locally func- 
tions of time and the independent spatial variables. 
For long-lived large-scale diabatic situations these coef- 
ficients vary slowly in space. Within mesoscale distur- 
bances they should vary in a pattern consistent with 
the disturbance. For example, the convection may 
moisten the upper atmosphere, giving a u and a v a pos- 
itive sign, while drying lower levels making the coef- 
ficients negative, Equations (12) and (13) allow a cli- 
matological estimate of the nonconservative property 
of the horizontal kinematic moisture fluxes to be easily 
determined from model calculations. A more complete 
solution procedure that makes no simplifying assump- 
tions will be presented elsewhere in an independent 
study. 

Substituting ( 12) and (13) into (10) and (11) gives 


-Jv g +fi> + u(A - a u q) = 0, ( 14) 

fu g -fu + v(A - a v q) = 0. (15) 


Note that A , a u q , and arf couple the two algebraic 
equations together. Dividing ( 14) and ( 15) by / and 
solving the coupled system for u and v yields 


_ Ug + MyVg 

“ 1 + M U M V ’ 


(16) 


Vg Ug 

1 + M U M V ‘ 


(17) 


Here the subscript m denotes the moist solution. Note 
that u m = u g and v m = v g when M U = (A- a u )f~ l and 
M v = (A — a v )f~ x are set to zero. Here A is negative 
for condensation, so a positive value of a u and a v would 
enhance the magnitude of M u and M v . For some cir- 
cumstances — for example, for large hourly condensa- 
tion rates in a hurricane — the magnitude of the non- 
dimensional number Ri = A f~ ] can exceed unity. A 
condensation rateof£ -0.8 X 10 _7 g g _I s _I yields 
R\ = -0.408, provided /= 10~ 4 and an average q 
satisfies q - 2 g kg -1 . Magnitudes for R x of 0.2 and 
less, however, are probably more common for large 
areas. 

For tropical flow the moist wind is computed using 
(6), (7), (12), and (13). This yields 



provided neither a u — A nor a v - A is equal to zero. 
The perturbation in the horizontal velocity is found 



July 1993 


RAYMOND 


1996 


using the pressure gradient computed from the diabatic 
heating. If the total derivative term; that is, in (6) and 
(7) is small, then the perturbation flow induced by the 
diabatic process is antitriptic; that is, there is a near 
balance between the Rayleigh damping term and the 
pressure gradient term. An example of antitriptic flow 
is found in Raymond ( 1986), illustrating topograph- 
ically induced mesoscale circulations in barotropic flow 
around an isolated mesoscale obstacle. Similar circu- 
lations containing positive and negative couplet pairs 
of convergence-divergence and cyclonic-anticyclonic 
flows are produced by this antitriptic mechanism. 

For extratropical flow an examination of ( 16) and 
(17) expanded in powers of provides a clearer pic- 
ture. Assuming that u m = u L + u^ n and v m - v L 
+ tJpcrt, where u L and v L are the horizontal wind com- 
ponents in the large-scale flow, and assigning all terms 
involving Q to the perturbation, yields for the pertur- 
bation and large-scale flow the expressions 


scale flow by factoring out V L but replacing all occur- 
rences of V L in the denominator by V g yields 



1/2 

, (22b) 


where V g = ( u ] + v\) ul . Substituting from (20) and 
(21 ) into (22b) gives 


V = V 

v m v g 


j _ ^2 + 2j?|(q u + q„) 


2Ri(ct u Ug 4* q u U^) UgVg(oc v ff K ) 


fv\ fv ? 


2RiUgVg{o[y q„)(q C i 

pv\ 


+ 


1/2 


Mperl = R\V g ~ R]u g + — 


R)(c*u + a v )u g 


U L = U g 


f 

Ufj , >Vg ^ (X u (X v Ug 


+ •••, (20a) 


f P 


+ ■ • (20b) 


^pert R\U g R i Vg + 


R i ( "F ) Vg 
— 


(21a) 


.. — «. i <*uUg ot u a v v g 

V L- Vg + — Jj— + • • *. (21b) 


Note that the large-scale flow is nearly balanced if 
\ct v \lf< 1 and | a u \/f < 1. Also the series in (20b) 
and (21b) converge rapidly, provided these ratios are 
much less than 1 . When this is not the case, this series 
expansion should not be used. 

In the Northern Hemisphere, given R } < 0 (con- 
densation) with u g > 0 and v g > 0, (20a) and (21a) 
imply that the magnitude of a westerly flow is reduced, 
since u g /u m > 1, while the magnitude of a southerly 
flow is enhanced, since v g /v m < 1. Repeating this 
thought experiment with an easterly flow leads to the 
reverse conclusions. In the Southern Hemisphere, u g 
> 0 implies that a northerly flow is enhanced because 
the coriolis parameter changes sign making R { > 0. As 
will be shown in the hurricane example given below, 
the diabatic perturbation terms [(20a) and (21a)] add 
significantly to the local convergence. 

The diabatic contribution alters the speed V m 
= (Um + v 2 m ) l/2 . Writing the definition of speed in 
terms of the large-scale wind and the diabatic pertur- 
bation gives 

Vm = W\ + 2 Hi I/pert + 2V L V pcn 

+ (Wper,) 2 + (Upcrt) 2 ] 172 - (22a) 
Here V L = (ul + vl) l/2 . Approximating the large- 


(22c) 

Note that V m = V g ( 1 — /?i) 1/2 when a„ = q„ = 0. All 
terms in the brackets except for the first two tend to 
nullify each other or make small contributions pro- 
vided \a v \/f<£ 1 and | a u \/f<£ 1. Thus, (22c) implies 
that the speed associated with moist cyclonic activity 
is less than the geostrophic value. 

For comparison purposes, note that the relationship 
between the geostrophic and gradient wind speed V v 
is 



where R is the radius of curvature (Holton 1979). For 
cyclonic situations ( fR > 0) the ratio V & /V g is less 
than 1. Thus, the speed is reduced from geostrophic 
values by the gradient wind primarily because curvature 
effects are included. In regions having condensation, 
the rate needed to sustain a specific curvature can be 
determined approximately by equating (22c) and (23), 
and solving for Q . Approximating the solution by the 
first two terms in (22c) gives 



As an example, if /= 10 -4 , V/fR = 0.2, and q = 2 
gkg _1 ,then Q= -1.1055 X 10~ 7 gg -1 s _1 . Note that 
Q is proportional to fq . 


3. An example 

Idealized flow in a hurricane is used to illustrate the 
diabatic influence on the wind field. An idealized axi- 
symmetric vortex flow was selected because of its in- 
teresting theoretical properties and because of the ini- 
tialization difficulties encountered when inserting a 
diabatic vortex flow in a numerical model. Following 
the procedure presented in section 2, but for cylindrical 
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coordinates and for axisymmetric flow, allows the moist 
radial wind component u m to be written in terms of 
the tangential velocity v according to 

u m = +jj . (25) 

It is assumed that the contribution from a v f~ l in M v 
is small. Near the vortex core v/fr P 1; thus, u m sim- 
plifies to 

u m — frM vy (26) 

while far from the core v/fr < 1, so 

u m = vM v . (27) 

The divergence is 



In the absence of condensation or evaporation (Q 
- 0) our solutions satisfy u = 0 unless a v + 0. Thus, 
radial velocities associated with condensation bands in 
hurricanes should be localized and exist in addition to 
a dynamically induced large-scale convergence asso- 
ciated with a v + 0. These bands in the radial velocity 
are seen in both numerical model simulations and in 
observations of hurricanes (Willoughby 1979), 

With a balance between the centrifugal and pressure 
gradient forces that characterizes intense atmospheric 
vortices, cyclostrophic flow leads to a tangential velocity 
defined by 



Therefore, solving for the moist tangential velocity v m 
= u, in the cylindrical moisture flux equations, using 
a procedure analogous to that developed in earlier sec- 
tions, gives 

v m “1“ ffVm Jru m Kf u — v c . (31) 


For v/fr P I, the divergence [(28)] simplifies to 


dr r 


« 2/A/, + 


rfdM v 
dr 9 


(29a) 


while when v/fr 1 the expression is 



In (29a) the sign depends upon M v and its radial vari- 
ation. In (29b) the divergence depends upon the vor- 
ticity (dv/dr + v/r) and the radial variation in M v . In 
both (29a) and (29b) M v < 0 contributes toward the 
convergence, provided the magnitude of M v is increas- 
ing radially outward and v > 0. Positive vorticity, av- 
eraged over the hurricane inner core, of 10 -2 s' 1 is 
common. Thus, as indicated by (29), substantial con- 
vergence is generated by diabatic processes when M v 
is sufficiently large in magnitude and of negative sign. 
Negative vorticity with M v < 0 and decreasing outward 
implies divergence. Evaporation near the cloud tops 
also encourages divergence. 

To express the radial velocity as a function of dis- 
tance from the vortex center, it is necessary to know 
M vy or Q when a/* 1 is small, and its radial distri- 
bution. Near the vortex center a constant value for M v 
would imply, following (27), a flow that is proportional 
to r (solid rotation), while if M v were linear in r, then 
u would be proportional to r 2 . Far from the vortex 
center the radial component [(28)] is proportional to 
the tangential component v and to M v . In the outer 
reaches of a hurricane the radial distribution of v de- 
creases approximately as r~ U2 (Riehl 1963). Thus, the 
radial variation in u m is considerable but depends a 
great deal upon the radial variation of the net conden- 
sation rate Q. 


Ignoring the dependence of u m on v m yields a moisture- 
modified gradient wind expression 


fr fr 
v m = — — ± — 

m 2 2 


J _ 4i/ m Af M + 4u 2 


11/2 


Jr (fr) 2 \ 


(32) 


Retaining the dependence of v m on u m developed above 
[(25)] requires solving the cubic equation 

vi + 2 :frv 2 m + [( 1 + M v M u (fr) 2 - v?]u m = frv] 

(33) 


for v m . Note that (33) remains cubic even when M v 
= M u = 0, implying than an extraneous real root 2 has 
been introduced by the algebra. 

Results are presented in Fig. 1 to illustrate the dif- 
ferences between (32) and (33). The axisymmetric 
tangential flow is taken from the analytical expression 
given in (6)ofFiorino and Elsberry ( 1989), with max- 
imum velocity F max = 35 m s _l . In their formula for 
the tangential wind the radius of maximum tangential 
winds is set to r max = 100 km and the parameter b 
= 0.96. For simplicity, in our experiments this radial 
profile is assumed to be in cyclostrophic balance, and 
is illustrated by the thick line in Fig. 1. Adding the 
influence of the Coriolis parameter (gradient wind) for 
a dry case [u m = 0, (32)] produces reductions in the 
tangential velocity of approximately 4 ms -1 at r mju , 
with larger values outside the radius of maximum 
winds. Willoughby ( 1990) indicates that the gradient 
balance condition closely approximations the axisym- 
metric tangential winds observed in hurricanes. How- 
ever, it should be remembered that measured pressure 


2 For the used, ( 33) has three real roots, two of which are negative. 
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Fig.-I. Idealized axisymmetric hurricane cyclostrophic gradient 
and moist tangential velocities (m s’ 1 ) as a function of radius. Speed 
and radial component also illustrated for the moist wind 


gradients already contain diabatic contributions. This 
could explain the discrepancy between these observa- 
tional findings (Willoughby 1990) and the theoretical 
and numerical findings that intense vortices react 
to specified heat sources (Hank and Schubert 
1986; Delden 1989). Note that our moist wind solu- 
tion, as illustrated by the divergence in (29), is also 
dependent upon the intensity of the vortex. 

To solve (33) for v m , the nondimensional values of 
M v and M u (Fig. 2) are assigned to be proportional to 
v c so that each is -2 at r mai . This large magnitude for 
u and M u is used to dramatize the response of the 
drag terms. The roots of (33) give one real positive 
solution for v m that can be used to evaluate (25 ), yield- 
ing a radial inflow of approximately 12 m s -1 . As ex- 
pected, the diabatic influence reduces the speed and 
shifts the radius of maximum tangential velocity from 
100 to 80 km, while generating the maximum mag- 
nitude of the radial velocity at r = 140 km. The max- 
imum reduction in the tangential velocity component 
is about 8 m s’ 1 , but the overall reduction in the speed 
is less since a strong radial inflow is present. Considering 
that our estimates of M v and M u are inflated in mag- 


nitude and in areal coverage, we can probably infer 
that the tangential velocities are not greatly altered by 
the diabatic drag term. 

Magnitudes of AC and M u dictate the radius of max- 
imum heating (RMH). Figure 2 shows a curve for M 
= M V = M u where the RMH is at 100 km, and another 
M that is positive between 0 and 50 km, illustrating 
possible modifications (exaggerated) because of inner 
core evaporation. These A/’s are used in the calculation 
of u m to show the importance of the divergence as a 
function of the RMH, as illustrated in Fig. 3. In every 
case the balanced solutions show that the greatest con- 
vergence is located interior to the circle with the pre- 
scribed RMH. Whether this would result in a change 
in the location of the RMH depends on the availability 
of moisture and other factors. However, in these cal- 
culations the radius of the maximum cyclostrophic 
wind was held fixed at 100 km. Otherwise, little change 
occurs in the magnitude or location of the maximum 
value of v m as a function of RMH. 

Observations or results from model simulations of 
intense vortices are needed to complete the evaluation 
of the moisture-modified horizontal winds and the 
evaluation of the diabatic and total derivative terms. 
The contribution from the total derivative would be 
expected to vary spatially, but the pattern should be 
consistent between similar types of situations and 
should change in time somewhat slowly. 

4. Applications 

Cumulus initialization and parameterization of long- 
lived convection in the extratropics are areas that could 
benefit from the knowledge contained in (16), (17), 
(18), and (19). The following discussion highlights 
some simple tests to investigate this potential in cu- 
mulus parameterization. These tests are presented to 
illustrate possible applications and do not represent a 
completely definitive study of the topic. 

Cumulus parameterization schemes — for example, 
Kuo (1965, 1974) — estimate subgrid-scale moist con- 
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F,G ' 2 ' Rad ' al Profiles of M illustrating (exaggerated) radial 
distribution of the condensation or evaporation rate. 
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Fig. 3. Radial distribution of the moist divergence (s _l ) 
as a function of M. 
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vective activity using the instantaneous state of the grid- 
scale variables and moisture or cloud properties. This 
information is then used to modify the large-scale tem- 
perature and moisture fields. Because cumulus clouds 
contain a large vertical component of motion, relative 
to their environment, they also are responsible for sig- 
nificant vertical momentum transport. Several schemes 
based on cumulus friction concepts have been intro- 
duced to adjust the grid-scale momentum (Ooyama 
1971; Shapiro and Stevens 1973; Schneider and Lind- 
zen 1967; Shapiro and Stevens 1980; Zhang and Cho 
1991). These approaches introduce a cloud velocity 
that-requires an additional determination, usually from 
a cloud model. In contrast, the terms in (20a) and 
(21a) that represent approximate changes in the large- 
scale wind introduce no cloud-scale dependent vari- 
ables but require calculation of the geostrophic wind. 

To test the momentum adjustment suggested by the 
perturbation terms in (20a) and (21a), they were added 
to the Kuo cumulus parameterization scheme in the 
local [Cooperative Institute for Meteorological Satellite 
Studies (CIMSS)] version of the Australian Bureau of 
Meteorology Research Centre’s (BMRC) finite-area 
model (Leslie et al. 1985). Thus, w ncw = w 0 id + Wpen 


where (20a) defines w^rt; similarly, u ncw = t?ow + v pcn 
where (21a) defines Our configuration of this 
semi-implicit regional model has 15 sigma levels. A 
12-h control forecast, 0000-1200 UTC 22 April 1988, 
with 1 50-km horizontal resolution and 600-s time steps 
was performed using the standard Kuo cumulus pa- 
rameterization scheme ( 1965) without any momentum 
adjustment. In the test forecast the first two terms in 
(20a) and (21a) are used to modify the horizontal wind 
components u and v, with the added requirement that 
| Ri | < 0.5. (Smaller magnitudes for may be more 
acceptable.) This large magnitude for was used to 
highlight the differences between the two schemes in a 
short forecast. Even though the convection may vary 
slightly in location and intensity between model time 
steps, it is still assumed that the convection is long 
lived. Spatially isolated convection will not result in 
introducing numerical noise because such noise is re- 
moved efficiently by the high-order low-pass tangential 
or sine filters (Raymond 1988; Raymond and Garder 
1991). In Fig. 4 the differences between the forecasts 
for the 500-mb mixing ratios (g kg" 1 ) are superimposed 
upon the infrared satellite image. Note that small in- 
creases in the mixing ratio occur over cloudy regions, 



Fig . 4. Difference ( x 1 0 ) in the mixing ratio (g kg ~ 1 ) for a 1 2-h forecast (0000- 1 200 UTC 22 April 1 988 ) 
with and without the perturbation wind [Eqs. (20a) and (21a)] application in the cumulus cloud param- 
eterization scheme. Differences are superimposed upon the VAS infrared image. 
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but decreases occur over regions associated with de- 
scending vertical motion. In this figure the band of 
clpuds to the north and west is associated with a cold 
front, while the band over the Qulf coast is high clouds 
streaming in over Mexico from the Pacific Ocean. 
(These bands are separated by a strong wind jet.) Our 
tests show that the. first term in (20a) and (21a) is 
responsible for most of enhanced vertical circulation. 

The convective rainfall from the cumulus scheme 
with momentum adjustment was slightly less than that 
produced by the traditional Kuo procedure because 
the restriction | R, | < 0.5 limited the amount of con- 
densation. This made the parameterization more re- 
alistic since too much condensation was occurring too 
rapidly in some localities with the standard scheme. It 
is known — for example, Krishnamurti et al. ( 1980) — 
that the original Kuo (1965) scheme gives excessive 
moistening of the atmosphere. Sundqvist et al. ( 1989) 
introduced a special factor to govern the vertical dis- 
tribution of convective heating to help correct that 
problem. Nevertheless, a maximum increase in the 
nonconvective rainfall of 4 mm in 12 h was also ob- 
served. This increase (Fig. 5) occurred because of the 
enhanced vertical circulation that resulted from the 


application of the momentum adjustment [(20a) 
and (21a)]. 

Initialized National Meteorological Center global 
analyses valid over North America between 1200 UTC 
9 March and 1200 UTC 12 March 1988 were used as 
initial and boundary conditions in a 72-h forecast made 
by the CIMSS regional model using 190.5-km resolu- 
tion on a 41 X 37 horizontal grid with a 720-s time 
step. The control forecast uses no momentum adjust- 
ment. In the test case the momentum adjustment uses 
the first two terms given in (20a) and (21a) with the 
restriction that | R, \ < 0.25. Also, the sixth-order im- 
plicit tangent filter (Raymond 1988) is applied with a 
filter factor of either < = 0.05 or « = 0.64. The latter 
value gives moderate to heavy filtering. The u velocity 
component rms differences between verifying initial- 
ized analysis and forecast are shown as a function of 
time in Figs. 6a-c. The momentum adjustment at 
sigma level 3 (approximately 230 mb) results in a small 
improvement in the statistics except at 24 h. Increasing 
the filter factor from 0.05 to 0.64 has a larger impact 
than the adjustment, primarily because all grid points 
are involved. In contrast, the maximum number of 
occurrences for cumulus parameterization is no more 



Fig. 5. Same as Fig. 4 but illustrating the 12-h difference in the total rainfall (mm) 
between the two 12-h forecasts. 
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Fig. 6. Root-mean-square differences are illustrated for the u ve- 
locity wind component at (a) sigma level 3, (b) sigma level 6, and 
(c) sigma level 10. “Modified” means the momentum adjustment is 
used. These statistics are computed with an epsilon filter factors (ep) 
of 0.05 and 0.64 in the implicit tangent filter. 


Time 


than 527 vertical profiles out of about 1 500 possible 
locations (see Fig. 8). At sigma level 6 (Fig. 6b, ap- 
proximately 440 mb) the major improvement occurs 
late in the forecast. Lower levels (approximately 720 
mb) respond best to the momentum adjustment, as 
illustrated in Fig. 6c. In all cases, heavy filtering always 
reduces the contribution from the adjustment. The v 
velocity wind component rms differences are shown in 
Fig. 7. Again, sigma level 10 shows the most improve- 
ment at 60 and 72 h. Upper-level contributions average 
out to near neutral or slightly on the negative side. This 
may reflect shortcomings in the cumulus parameter- 



Time 

Fig . 7. Same as Fig. 6 but for the t? velocity wind component 
and with the smaller filter factor. 


ization or in the momentum adjustment process or in 
both. The lack of entrainment or vertical mixing by 
turbulence in the cumulus parameterization is clearly 
a deficiency. How case dependent the adjustment pro- 
cedure is remains unknown, as do estimates for a u 
and a v . 

The number of active cumulus parameterizations 
versus time step is shown for a 36-h period in Fig. 8. 
It shows that cumulus convective activity is slightly 
enhanced early in the forecast and experiences less 2 At 
and other high-frequency oscillations with the mo- 
mentum adjustment. 



Time Step 


Fig. 8. Number of grid points experiencing cumulus 
parameterization as a function of the time step. 
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5. Summary 

A diabatically induced term is present in the equa- 
tions describing the moisture flux and heat flux. This 
term has the form of damping (forcing) for conden- 
sation (evaporation). Diabatic processes also influence 
the horizontal flux through the pressure gradient term. 
In this study the influence of the damping term is eval- 
uated, but the characteristics associated with diabatic 
changes to the geostrophic wind are not analyzed. 

This study finds that mild reductions in the large- 
scale wind speed caused by cumulus friction can con- 
tribute to the near conservation of the large-scale hor- 
izontal components of the moisture flux. This is the 
consequence of an establishment of a near balance be- 
tween the pressure gradient, Coriolis, and Rayleigh 
damping terms in the moisture flux equations. This 
balance is identified by scale analysis. This outcome is 
independent of the mechanism responsible for the re- 
duction in wind speed. 

Estimates of diabatically related changes to the hor- 
izontal wind are easily obtained from the balanced re- 
lationship. For realistic condensation rates, these wind 
estimates are consistent with the gradient wind ap- 
proximation. When these estimates are used in a simple 
hurricane model they predict a radial inflow with con- 
vergence, provided the rotation is cyclonic. In contrast, 
divergence is associated with anticyclonic rotation and 
evaporation. When the predicted wind perturbations 
are used as a momentum adjustment in a Kuo cumulus 
parameterization scheme they improve the 72-h fore- 
cast rms u velocity component error statistics and en- 
hance mesoscale vertical circulations. In the tropics 
any movement toward conservation of moisture flux 
enhances an antitriptic relationship, cementing a near 
balance between a Rayleigh damping term and the 
pressure gradient term. Note that it is necessary to have 
the condensation before the wind-related processes can 
be activated. In this way this process is similar to con- 
ditional instability of the second kind (CISK). 

Our equations (20a) and (21a) are approximations 
of the influence from one diabatic source on the wind. 

A more rigorous study using the vorticity and diver- 
gence concept will be described elsewhere. The terms 
in the flux equations are also being analyzed using me- 
soscale model calculations. This is the subject of con- 
tinuing study. 
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APPENDIX 


The Sensible Heat Flux Equations 


The equations describing the temporal and spatial 
distribution of the horizontal components of the sen- 
sible heat flux can be constructed from ( 1 ), (2), and 
(4). They are 


d(8u) 

dt 


8 dp 

-~ + BJv- B8u, 
p dx 


d(9v) 

dt 


8 dp 
p dy 


Bfu - B9v. 


(Al) 

(A2) 


In these equations the damping coefficient B = LQ/ 
( Tc p )\ here, T is temperature, c p is the specific heat at 
constant pressure, and L is the latent heat of vaporiza- 
tion. For condensation processes B is negative but 
smaller in magnitude than the coefficient A in (8) 
and (9). 
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1 . INTRODUCTION 

In this study we examine the initialization and 
assimilation of cloud and rainwater quantities in a mesoscale 
regional model. In our study forecasts of explicit cloud and 
rainwater are made using conservation equations. The 
physical processes include condensation, evaporation, 
autoconversion, accretion and the removal of rainwater by 
fallout. These physical processes, some of which are 
parameterized, represent source and sink terms in the 
consei vation equations. Flic question of how to initialize the 
explicit liquid water calculations in numerical models and 
how to retain information about precipitation processes 
during the 4-0 assimilation cycle are important issues that 
will be addressed. 

With increasing computer speed there has been a slow 
but study enhancement in model horizontal grid resolution. 
Current N\V P models are approaching the spatial resolution 
where limc-dcpcndcnl cloud variables are physically relevant. 
Tn our pilot study we examine the feasibility of initializing and 
assimilating cloud and rainwater using explicit cloud 
conservation equations within a regional mesoscale model. 

I he knowledge of cloud and rainwater fields as a function of 
time and space provides information, such as model heating 
rates and vertical distributions, that are more difficult to 
obtain from conventional cloud parameterization approaches 
where cloud properties are not retained from one time step to 
the next. Global models will soon routinely solve explicit 
cloud and rainwater conservation equations, but limited area 
mesoscale models can be utilized to examine this topic now. 

Cloud and rainwater are the cumulative products of 
the transport of atmospheric moisture, both vapor and liquid, 
plus the physical processes included in cloud development 
and decay. By their very nature cumulative processes are 
difficult to initiate in numerical models, and it is common for 
precipitation to be absent during the first few hours of a 
numerical forecast. One reason for this deficiency is the 
resolution of the data. Conventional surface and upper air 
reports are known not to capture the mesoscale circulations 
that have the enhanced convergence necessary to support 
mesoscale precipitation processes. The missing mesoscale 
circulations also influence the moisture and thermal fields 
(Turpeinen 1990). Consequently, in order to initiate 
precipitation and assimilate cloud parameters it is necessai^ to 
know additional information beyond the conventional 
synoptic data base. This might include accurate estimates of 
the vertical heating rates and the low level convergence and/or 
the high level divergence that produces the cloud-sustaining 


vertical motion field. Also, model generated inaccuracies or 
shortcomings, c.g., phase errors and inadequate physics, 
must be accommodated in some manner. 

When the heating rales are known or can be 
approximated, then cither diabatic or cumulus initialization is 
possible (Errico and Rasch 1988; Donner 19S8). These 
initialization procedures help establish circulation patterns 
consistent with latent heating processes that are otherwise not 
resolved in synoptic scale data. For a review of recent 
studies using diabatic initialization schemes, see Table I in 
Turpeinen et al. (1990). Without the heating rate 
information each forecast must go through a short time period 
where the precipitation processes must ’spin-up 1 before 
rainfall can be produced. In regional mesoscale models this 
spin-up usually occurs during the first 2 to 6 hours, 
Diabatic initialization should reduce the degree and duration 
of the spin-up. In contrast, normal mode or dynamic 
initialization procedures, designed to help with the removal of 
gravity waves during initialization, only provide slight or no 
assistance in speeding up precipitation processes (Lejcnas 
1980). 

Cloudwaier and rainwater data arc generally not 
available. Rain gauge reports imply information about the 
total amount of liquid water reaching the surface, but no 
knowledge about the instantaneous vertical distribution is 
recoverable. Most microwave satellite measurements of 
instantaneous cloud quantities only resolve the vertically 
integrated cloud and/or integrated rainwater. Thus it is 
necessary to make some assumption about the vertical 
distribution of the heating rate. New physical retrieval 
methods that combine observed microwave satellite radiances 
with known meteorological fields and/or physical constraints 
should provide some insight into the vertical partitioning of 
the liquid water. These techniques are just beginning to be 
examined. 

The use of diabatic initialization requires an answer to 
the question of how the latent heating is going to be 
distributed in the vertical. Some studies have used 
representative profiles (Fiorino and Warner 1981); others 
assume a general parabolic profile (Tarbcll et al 1981; Salmon 
and Warner 1986, etc.) while some have tried to use the 
forecast model itself in some fashion (Wergen 1988: Danard 
1985: Turpeinen et al. 1990). Donner (1988) used a 
variational formulation to optimize the vertical distribution 
with respect to the methodology of the Kuo Cumulus 
parameterization. In all but t lie Donner study it has been 


customary to use either the diabatic omega equation or a 
diabalic form oT the normal mode initialization scheme. 
Earlier pioneering work by Krislinanuirli et al. (1984) 
suggested that a rcanalysis of the moisture and wind field is 
also needed. 

In this study we will incorporate in the initialization 
the heating rates, cloudwater, and rainwater predicted by a 
forecast beginning 12 hours before die assimilation process. 
This information will be combined with new information 
from observations, and heating rates will be utilized in a 
diabatic form of the vertical mode initialization scheme of 
Bourke and McGregor (1983). Theoretically, if data is 
always available then cloud and rainwater could be treated 
much like other meteorological variables, i.e., updated and 
modified during the 4D assimilation process. Caution is 
however required when making large changes to existing 
model cloud fields, or to the dynamics that support the 
precipitation process. Testing and experience with explicit 
cloud initialization will dictate what the correct procedure 
should be. 

Results from only a few experiments are presented 
lierc. An evaluation of the contribution of the satellite data 
will be given at the conference. 

2. MODELING CONSIDERATIONS 

A modified version of the Bureau of Meteorology 
Research Centre (BMRC) operational incremental O/I 
analysis (Mills and Seaman 1990), vertical mode initialization 
or VMI (Bourke and McGregor 1983), and semi-implicit 
forecast model (Leslie, et al. 1985) are used in this 
investigation. The incremental analysis is set to preserve the 
irrotational wind Held while the rotational component of the 
differences between data and model guess wind and other 
fields are evaluated by the O/I. In our initial tests the 
horizontal resolution is 120 km, while the vertical sigma 
coordinate is subdivided into 15 levels. We selected the low 
resolution 120 km grid spacing on a 35x35 horizontal grid 
configuration to specifically examine problems that are 
relevant or commonly encountered in today's centers for 
numerical weather prediction. 

In our present study the Kuo cumulus 
parameterization is turned off and the conventional large scale 
precipitation procedure is replaced with explicit conservation 
equations for cloudwater and rainwater. The explicit cloud 
physics contains processes modeled using a modified Kessler 
type scheme (Kessler 1974), while the condensation and fall 
velocity calculations are similar to that described in Amhes, et 
al. ( 1 987). The 'fallout' term is calculated by a semi- 
Lagrangian method to allow for the 10 minute time steps used 
1,1 ,hc semi-implicit model calculations. Ice phase cloud 
particles have not yet been included in the model. 

Part of the goal of the current study is to modify the 
model cloud and rainwater fields using SSM/I retrievals of 
integrated total liquid water or vertical distributions of water 
These observed quantities will be used to perturb existing 
model cloud and rainwater fields. Heating rates, lagged by 
one time step (At=!() min) or based on the total liquid water 
rom the I2li pre-forecast, are used in the initialization, 
any different types of model experiments are to be 


conducted, including those related to the model's response to 
diabatic and initial information. 

Because ol horizontal scale considerations, our 
explicit cloud calculations do not mode! individual clouds but 
do approximate average values for a system or complex of 
clouds. This was the motivation behind the modifications 
made to the Kessler (1974) cloud physics formulation. It is 
known that the lifetime of individual clouds is less than one 
hour while a cloud complex exists for many hours with 
remnants lasting for days. Thus Kessler's accretion 
(collection) process was modified to slow greatly the 
conversion from cloud to rainwater when the cloud water C is 
below a preset critical value C cri| . Also incorporated is a 
vertical dependence on temperature that follows the scheme 
suggested by Schlesinger (1990). This latter change helps 
retain and enhance anvil or cloud top features. In mixing 

ratio units (gg '), the contribution to the rale of change from 
accretion is 

ACC=0.8SC,C ; C R ns7s (I) 

where C,=min[ 1.0, 1.0 - (C crjl - C)/C ctU | , 

C 2 =min[ 1.0, (T-233)/40 | , 

C 2 -max( 1 .0, C 2 | . 

I fere R is rainwater in mixing ratio units, T is the temperature 
in degrees Kelvin, and C„. =1 gkg'. Initial experiments 
give promising results, but additional tests and verifications 
will be made to further evaluate computed cloud water 
quantities. These will include comparisons between model 
cloud water fields and SSM/I retrievals of cloud water 
distributions. 

To facilitate the large time step size used in the semi- 
implicit forecast model the fallout term in sigma coordinates is 
rewritten as two terms, one of which represents vertical 
adveclion, 

g0(pv,R)/i)o s AO( P *R)/0 ct + (p* R g/p * )D( p v t )/Da . (2) 

Merc p is the air density, V[ (ms 1 ) is the terminal velocity of 
raindrops, p* is the surface pressure, g is the acceleration of 

gravity and A (s ') is equivalent to a velocity in the vertical o 
coordinate system, that is 

A= gpVj/p . (3) 

The time tendency for A using (2) now can be iound by 
semi-Lagrangian procedures. 

3. DISCUSSION OK RESULTS 

Super-imposed over the GOLS 7 IR image in big. l 
are contours of cloudwater (solid) and rainwater (dashed) for 
sigma level 4 (approx. 300 mb) from a I2li forecast. The 
contour interval is 0.1 gkg The initial field for this ERICA 
(Experiment on Rapidly Intensifying Cyclones over the 
Atlantic) IOP-5 case is obtained from NMC's initialized 
global analysis. For the most part, the forecasted cloud 
fields assume the horizontal comma-shape shown in the 



Fig. 1. Contours (0.1 gkg' 1 ) of cloudwater (solid) and 
rainwater (dashed) from a 12h forecast (0000 to 1200UTC, 
20 Jan 1989) are superimposed over the GOES IR image. 

image. In this forecast no diabatic initialization was used. 
When the modifications to the accretion term proposed in Ft]. 
(1) are removed, the amount of cloud water is reduced to just 
one short contour (not shown) while the rainwater is 
enhanced slightly beyond that shown in Fig. I. Cloudwater 
amounts at lower sigma levels also are much reduced when 
the original Kessler formulation is used. With cither method 
the total rainfall is nearly identical in magnitude and areal 
coverage. The maximum rainfall over 12 hours was 33 mm 
for the modified approach and 34 mm for the original Kessler 
formulation. 

In Fig. 2 the area averaged precipitation rate at each 
lime step is plotted versus the time step for forecasts between 
1200UTC. 19 Jan and 1200UTC, 20 Jan 1989. In the two 
forecasts without explicit clouds (labeled ’’large scale”), all 
excess water vapor exceeding saturation falls out 
immediately. This approach is known to produce rainfall 
rapidly. In contrast, note that the precipitation in the explicit 
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Fig. 2. Area average precipitation rate (cm) is plotted 
versus the lime step for 4 forecasts. Forecasts with and 
without explicit clouds illustrate the differences in the 
precipitation spin-up time. 


cloud calculations is delayed. This is because the 
autoconversion threshold, used to model the water holding 
capacity of clouds, is not triggered until the cloudwater 
exceeds 0.5 gkg' 1 . 

After 12h or 72 time steps, the forecast fields 
generated by the explicit cloud calculations are assimilated to 
begin a second 12-hour forecast. Note in Fig. 2 that the 
precipitation rate remains unchanged following the 
assimilation step as illustrated by the continuous curve 
between steps 72 (square) and 73 (circle). However, the 
precipitation rale can incur irregular behavior if the irrotational 
wind or dynamics are significantly changed. The incremental 
analysis in our assimilation preserves the irrotational wind 
component from the previous forecast. But, our current 
VMI initialization can alter the convergence and/or 
divergence. Thus, some consideration is being given to 
installing an incremental initialization procedure. 

In Fig. 3 we expand the view of the second !2h 
period presented earlier in Fig. 2. Two of the curves from 
Fig. 2 are reproduced for comparison purposes. They 
include the cold started large scale calculations represented by 
the triangle curve, and the assimilated hot start shown with 
circles. Additionally, cold and hot started explicit cloud 
calculations beginning with the global initialized analysis are 
presented. The hot started forecasts utilize cloudwater and 
rainwater forecasted by the !2h pre- forecast. The heating 
rate at time step 72 is also used in the diabatic vertical mode 
initialization. 



Fig. 3. Same as Fig. 2 except for initialization. All 
forecasts begin at time step 72. 

It is clear that when the previous forecast model fields 
arc used, as in the assimilation process, there is dynamics to 
support the cloud and rainwater. Tlius the forecasted 
precipitation rale behaves properly (circle curve). When the 
supporting mesoscale dynamics and moisture supply arc 
missing, as in the global initialized analysis, then adding 
cloud and rainwater initially can result in less realistic 
behavior as illustrated by the heavy solid curve in Fig, 3. 
This curve illustrates our worst-case scenario. For this latter 
case, improvement is obtained when the heating rates are also 
used to force or nudge the dynamics (Wang and Warner 




1988) toward the desired goal (not shown). We discovered 
that enhanced development can become too strong in a 
reduced area when only the heating rates from time step 72 in 
the 12h pre-forecast are used. Approximating the heating 
rates from the total eloudwater and rainwater Helds gives 
better results. 

Many additional experiments are required before the 
best strategy or strategies are identified. But, the assimilation 
approach using incremental changes clearly is indicated to be 
a powerful technique that could solve the spin-up problem. 
However, remotely sensed eloudwater and rainwater data, 
and information about their vertical distribution, is still 
required if the forecast is to correctly simulate the 
atmosphere. 
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1. INTRODUCTION 

The retrieval of precipitation distributions by passive 
microwave techniques has been an area of study in satellite 
remote sensing for well over a decade. The recent initiative to 
estimate both the horizontal and vertical distributions of 
precipitating liquid water in the tropics (i.e. the Tropical 
Rainfall Measuring Mission or TRMM) has further stimulated 
interest in this area; see Simpson, et al. (1988). The proposed 
TRMM satellite will carry a suite of sensors including passive 
and active microwave instruments, as well as visible and 
infrared radiometers. A combination of microwave and 
visible/IR sensors has also been proposed for one of the Earth 
Observing System (EOS) platforms for the purpose of 
retrieving rainfall rates and other geophysical parameters; ref. 
Murphy (1987). Methods for determining precipitation 
distributions based upon data from several instruments having 
different spectral and resolution characteristics will therefore 
be required. 

Currently multispectral microwave observations from the 
Special Sensor Microwave/Imager (SSM/I), borne by the 
DMSP-F8, have provided researchers with the means for 
retrieving rainfall rates; see Hoi Unger (1991), Kummerow and 
Liberti (1990), Petty and Katsaros (1990). The SSM/I is a 
multichannel, dual-polarization, passive microwave radiometer 
with channels at 19.35, 22.235, 37, and 85.5 GHz. The 
diffraction limitation of the SSM/I antenna causes the spatial 
resolution of measurements to increase with frequency, such 
that the minimum footprint dimension decreases from 43 km at 
19.35 GHz to 1 3 km at 85.5 GHz. The DMSP also cajTies the 
Special Sensor Microwave/Temperature sounder with channels 
at 50.5, 53.2, 54.35, 54.9, 58.4, 58.825, and 59.4 GHz, 
and a footprint dimension of approximately 180 km at nadir 
view. In the present study, a physical retrieval method for 
estimating precipitation water distributions and other 
geophysical parameters based upon measurements from the 
DMSP-F8 SSM/I is developed. Three unique features of the 
retrieval method are (1) sensor antenna functions are explicitly 
included to accommodate varying channel resolution, (2) an 
embedded one-dimensional cloud model is utilized to generate 
vertical distributions of precipitating and nonprecipitating 
water, and (3) spectral solutions are sought for certain 
background parameters, such as humidity, which vary more 
slowly in the horizontal than the cloud and precipitation water 
contents. The general framework of the method will allow us 
to incorporate measurements from the SSM/T and 
geostationary infrared measurements, as well as information 
from conventional sources (e.g. radiosondes) or numerical 
forecast model Fields. 


Z RETRIEVAL METHOD 

The basis of the retrieval method is the minimum variance 
solution as described in Lorenc (1986). An error functional 
which expresses the deviation of the observed satellite antenna 
temperatures from model-derived values, plus an additional 
constraint which represents the deviation of retrieved 
geophysical parameters from a priori estimates, is derived. 


E = (TAobs - TA(P) m0 d) T 2 (TA obs - TA(P) mod ) 

+ 0* * Papriori)^ A ^ (P-Papriori) (1) 

Here, TAobs > s a vector of antenna temperature 
measurements, TA(P) mCK j are antenna temperature estimates 
based on a physical model, P is a vector of geophysical 
parameters to be retrieved, and Papriori are estimates of the 

parameters based upon a priori knowledge. £ and A are 
covariance matrices representing observational plus model 
errors (I), and errors in the a priori estimates (A). In this 
study all off-diagonal elements of these matrices are neglected. 

The antenna temperature model, TA(P) mo d» embodies 
the physical relationship between the geophysical parameters 
to be retrieved and the antenna temperatures measured by the 
radiometer. An element of TA(P) m od may be written 

TA(P) p = A c T TB(P) p + A x T TB(P) p ' + 5 p Tbb, (2) 

where Ac and A x are the co- and cross-polarized antenna 
patterns of the radiometer. TB(P) p and TB(P) p - are the 
modeled brightness temperatures in the polarization p and 
orthogonal polarization p' with respect to the plane of 

polarization of the measurement, and 8 p is the fraction of the 
radiometer feedhom pattern not subtended by the antenna. 
Tbb is the cosmic background radiance (2.7 K). Thus each 
antenna temperature is modeled as the convolution of the 
upwelling brightness temperature field by the antenna response 
pattern of the sensor. In this way, measurements at different 
channel frequencies or from different sensors which have 
different sampling / spatial resolution can be accommodated. 

The retrieval of preciptation liquid water contents and 
other geophysical parameters is accomplished by iteratively 
perturbing the geophysical parameters until the error functional 
Eq. (1) is minimized. Unlike Lorenc (1986), we minimize the 
error functional using the memoryless, quasi-Newton method 
described by Shanno (1978). 


3. SENSOR RESPONSE MODEL 

In this section, the components of the sensor antenna 
temperature model, Eq. (2), arc described. 

3.1 Model Grid 

Both the antenna response functions and brightness 
temperature fields are discretized on a model grid. Rectangular 
grid "boxes" are defined on the grid, such that the center of 
each box coincides with the earth location of an all-channel 
antenna temperature measurement from the SSM/I. The boxes 
are oriented along the SSM/I A-scans (which contain the all- 
channel measurements) at a regular spacing of 25 km. Since 



the cross-scan separation of succeeding A-scan lines varies 
along the scan line, the cross-scan dimension of the grid boxes 
also varies, with a maximum dimension of 25 km in the center 
of the swath, and decreasing towards swath edge. In the 
discretization of the antenna response functions and brightness 
temperature fields, the variation of the grid box dimension 
with scan position is taken into account 

3.2 Antenna Response Functions 

Measurements of the SSM/I antenna reponse functions 
and a detailed description of the SSM/I scan geometry were 
provided by Gene Poe (currently with Aerojet Electrosystems 
Co.). Using this information the integrated co- and cross- 
polarized antenna responses over each grid box for all SSM/I 
measurements are computed 

3.3 RadianvgJjansfer / Cloud Model 

The model for the average polarized brightness 
temperature, TB(P) p , upwelling from a grid box is described 
here. The earth's surface is characterized by a skin 
temperature and an- emissivity. The emissivity in each 
polarization is represented using the two-parameter formula of 
Grody (1988). The downwelling brightness at the top of the 
atmosphere is set equal to the cosmic background value (2.7 
K). 


The model atmosphere is divided into 50 mb layers, 
starting from an assumed top of the atmosphere at 100 mb. 
Atmospheric temperature falls off linearly with height at a 
fixed lapse rate of 5.4 K/km from a variable surface value. 
The surface relative humidity (and the associated water vapor 
density) is also allowed to vary in the retrieval. It is assumed 
that the water vapor density falls off exponentially with height, 
with a scale height of 2.0 km. The temperature lapse rate and 
vapor scale height reflect the tropical climatological values 
compiled by Jordan (1958). 

Clouds are assumed to occupy a variable cloud fraction 
within each grid box. Within the cloud fraction, the vertical 
profiles of precipitating liquid water and ice hydrometeors, as 
well as nonprecipitating cloud droplets and ice, are determined 
using a one-dimensional, steady-state cloud model. The 
updraft in the cloud model is an adaptation of the plume model 
described in Anthes (1977), which utilizes the water substance 
conservation equations of Kessler (1965). In light of the 
arguments by Schlesinger (1990), the rate of accretion of 
cloud water by rain drops is reduced by a factor of three from 
the original formulation in Kessler. For a given cloud 



percent; coastlines are dashed. 



Fig. 2. Retrieved liquid precipitation water contents in the 
vicinity of Emily at 10 UTC on September 21, 1987. Fields 
are shown at pressure levels of (a) 550 mb, (b) 700 mb, (c) 
850 mb, and (d) at the surface. Contours are drawn at 
increments of 0.2 g/m^, starting with a minimum contour of 
0.2 g/m^. Maximum values are in tenths of a g/m^; coastlines 
are dashed. 




environment, vertical development of the cloud is controlled 
by a single parameter, the horizontal mass entrainment rate. 
The entrainment rate is assumed to be inversely proportional _to 
the cloud updraft radius, as in Simpson and Wiggert (1969). 
The updraft radius is allowed to vary in the retrieval. The 
transition from liquid to ice phase water substance is linear in 
the model, beginning at the 273 K isotherm and ending at 253 
K, above which the cloud is completely glaciated. 
Precipitation fallout is assumed to occupy a separate downdraft 
region, which has the same horizontal cross-sectional area as 
the updraft. Therefore, within the area determined by the 
cloud fraction, the liquid water contents of the various cloud 
water constituents at each level are the averages of the updraft 
and downdraft values. 

Absorption coefficients for molecular oxygen and water 
vapor are computed using the formulae of Liebe (1985), and 
cloud absorption coefficients are computed based upon 
Rayleigh theory. Extinction coefficients, single-scatter 
albedoes, and asymmetry factors for liquid and ice 
hydrometeors are determined using fits to detailed calculations 
by Kummerow and Weinman (1988). 

Brightness temperatures upwelling from each grid box 
are computed using a multiple-scattering radiative transfer 
code. The radiative transfer solution is based upon 
Eddington's second approximation for a vertically-layered 
atmosphere. Separate radiances are computed for the clear and 
cloudy regions of the grid box, and then averaged with a 
weight equal to the fractional coverage of each region within 
the grid box. 


4. APPLICATION OF THE RETRIEVAL METHOD TO 

SSM/I DATA 

The physical retrieval method is applied to SSM/I 
observations of the incipient Hurricane Emily at 10 UTC on 
September 21, 1987, as it was developing over the Windward 
Islands off the northeast coast of South America. A total of 
eight geophysical parameters are optimized in the retrieval: 
four emissivity parameters (two for each polarization), surface 
air temperature, surface relative humidity, cloud fraction, and 
cloud updraft radius. Since only thirteen observations per grid 
box are available from the SSM/I, band-limited spectral 
solutions for the emissivity parameters, air temperature, and 
relative humidity are sought in order to greatly reduce the 
number of degrees of freedom in the problem. Only the DC 
components of the emissivity parameters (one for land and one 
for ocean regions) and surface air temperature are utilized to 
represent these parameters over the solution grid. A linear 
trend, the DC, and first wavenumber are used to represent 
horizontal variations of the surface relative humidity field. A 
separate cloud fraction and effective updraft radius for each 
grid box are also determined. Also, only the 37 and 85.5 GHz 
SSM/I data are utilized in the retrieval to be presented. 

After 15 iterations the error functional Eq. (1) without 
any additional constraints (second term equal to zero) is 
effectively minimized to produce the surface relative humidity 
and precipitation liquid water content fields in Figs. 1 and 2, 
respectively. The retrieved fields in this first test case appear 
reasonable, although the precipitation in the core region of 
Emily is extreme (equivalent surface precipitation water 
content greater than 1 g/m^ or 18 mm/hr covering a large 
area), and the vertical development of the associated cloud is 
limited. 


5 EXTENSION OF CURRENT WORK 
AND APPLICATIONS 

Work will continue to further improve the representation 
of clouds and their environments in the sensor response 
model. The relatively rapid convergence of the retrieval 


method (all fields had essentially converged in less than 10 
iterations) using only four channels of sensor data suggests 
that more degrees of freedom may be allowed in the retrieval 
solutions. Investigations by Short (1988) and others have 
indicated that the representation of clouds by a statistical 
distribution over each grid box rather than by a single vertical 
profile may be superior for remote sensing applications. 
Greater variability in the solution fields of emissivity, air 
temperature, and water vapor may also be allowed if additional 
sensor data and/or constraints are utilized in the retrieval. 
Future tests will include the lower-frequency SSM/I data as 
well as SSM/T and geostationary infrared measurements. 
Retrievals will be validated using volume-scan radar 
measurements and conventional weather data. 

One of the goals of the present investigation is to utilize 
retrievals of precipitation water contents as a proxy for latent 
heating in rapidly developing storms. Latent heating 
information may help to improve short-term numerical weather 
prediction model forecasts of storms, but forecast model fields 
may also be used to help improve retrievals through the 
constraint term in Eq. (1). The interplay between forecast 
models and retrievals is a subject of further study by the 
authors. 
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SSM/I Physical Retrieval of Precipitation Distributions 
and Numerical Model Forecast Applications 
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The Special Sensor Microwave/Imager (SSM/I) is a passive microwave radiometer with 
channels at 19.35, 22.235, 37, and 85.5 GHz. At the lower three frequencies, absorption 
and re-emission of microwaves by liquid precipitation produces a measurable signal which 
is strongly related to the vertically-integrated liquid water within the radiometer Field of 
view. At 85.5 GHz, liquid precipitation also has an absorption/emission signature, but the 
upwelling microwave radiances are strongly modulated by scattering off precipitation-sized 
ice particles which occur at greater altitudes. The combination of channels can therefore be 
utilized to infer approximate vertical distributions of precipitation water contents. 

A physical retrieval method based upon a minimum variance estimator is applied to the 
SSM/I measured radiances to deduce the vertical distributions of precipitation water 
contents in different weather systems. The retrieval method incorporates a radiative 
parameterization for computing upwelling microwave radiances at the SSM/I frequencies as 
functions of the non-homogeneous distributions of precipitating liquid and ice 
hydrometeors over a 25 km x 25 km area (the sampling resolution of the SSM/I). The 
radiative parameterization is based upon detailed calculations of the propagation of 
microwaves through simulated cloud fields. Both convective and stratiform cloud fields 
are generated using a numerical model. The spatially inhomogeneous response of the 
SSM/I instrument to the upwelling radiance field (an antenna diffraction effect) is also 
explicitly modeled in the physical retrieval method. 

Retrieved precipitation water distributions are compared to volume radar scans from 
Darwin, Australia and Kwajalein, Marshall Islands. 

Latent heating rates inferred from SSM/I retrieved precipitation water distributions can 
be assimilated into numerical weather prediction model forecasts. This kind of information 
is especially important in data-sparse regions (i.e., over the oceans) where the development 
of new storms may not be captured by the traditional observational network. Latent heating 
fields have been assimilated into CIMSS Subsynoptic Scale Model forecasts using a new 
method, and have shown significant impact on the "spin-up" and subsequent evolution of 
rapidly deepening storms. 

Latent heating information is used in the model's diabatic initialization and in the 
Newtonian nudging of the temperature field early in the model forecast. The nudging 
commences with the initial guess and is modified through time by the forecasted wind field 



and condensation rates. The Newtonian nudging of the model temperature field at locations 
experiencing latent heating represents an external forcing that enhances rapid deepening in 
model forecasts. Our findings show that diabatic initialization using two vertical modes 
gives a small response, while four vertical modes can produce large changes that on 
occasion generate model instability. Important factors influencing the outcome include the 
model grid size, mesoscale content of the initial model fields and the model physics. By 
introducing an estimated latent heating field into a numerical simulation of the ERICA 
IOP-5 storm, 24 hour forecasts of central pressure decreased as much as 20 mb. 
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1 . INTRODUCTION 

In the physical retrieval of precipitation 
amounts using the satellite-borne Special Sensor 
Microwave/Imager (SSM/I), an accurate radiative 
transfer model relating liquid and ice-phase 
precipitation water contents to upwelling, cloud top 
radiances is required. One such model that has been 
used to retrieve precipitation amounts from SSM/I (ref. 
Olson and Raymond, 1992) is based upon Eddington’s 
second approximation, which includes the effect of 
multiple scattering of microwave radiances off 
precipitating hydrometeors. In the present study radar- 
derived hydrometeor distributions will serve as input 
to this radiative transfer model, to determine whether 
or not consistency between the model-output radiances 
and SSM/I observations can be achieved, and under 
what modeling assumptions. This study should lead to 
an improved radiative model and more accurate SSM/I 
physical retrievals. 

2 . RADAR DATA AND PROCESSING 


Volume radar reflectivity data from the 5 -cm 
Darwin/TOGA installation recorded during the passage 



of a squall line at 18:20 UTC on February 9, 1988 are 
interpolated to a cartesian grid with dimensions 240 
km x 240 km x 18 km. The grid spacings are 2 km in 
the horizontal and 1.5 km in the vertical. The radar 
data in elevation/azimuth/range format are interpolated 
to grid locations based on the method of Mohr and 
Vaughan (1979). 

Within the convective leading edge of the 
squall line, the fractions of liquid and ice precipitation 
associated with a measured reflectivity are estimated 
using the temperature-dependent partitioning scheme 
of Moss and Johnson (1992). Within the trailing 
stratiform precipitation region, it is assumed that the 
radar reflectivity is due entirely to ice-phase 
hydrometeors above the freezing level, and below the 
6 C level the radar reflectivity is due entirely to liquid 
hydrometeors. A linear transition from ice to liquid 
between 0 C and 6 C is assumed. This stratiform 
precipitation partitioning scheme is based upon model 
simulations of anvil clouds using a one- dimensional 
version of the Ridout (1991) model. The vertical 
profile of environmental temperatures is obtained from 
the Darwin sounding at 22:00 UTC on February 9 
1988. 



Fig. 1. Fields of slant path- integrated liquid precipitation (a) and ice precipitation (b) derived from the 
Darwin/TOGA radar at 18:20 UTC on February 9, 1988. Coastlines are indicated by white lines, and the white 
circle indicates a data-void region in the radar volume scan. Each panel covers an area of 240 km x 240 km. 
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Fig. 2. Simulated brightness temperature fields at 37 GHz and 85.5 GHz are shown in panels (a) and (b) f 
respectively. Simulated antenna temperature fields at 37 GHz and 85.5 GHz are shown in panels (c) and (d), while 
SSM/I observed antenna temperatures at 37 GHz and 85.5 GHz are shown in panels (e) and (f). 






Integrated liquid and ice-phase precipitation 
fields derived from the Darwin volume scan are shown 
in Figs, la and lb, respectively. The convective 
leading edge of the squall line is seen crossing the 
Darwin radar site along the northern coast of Australia. 
The region of trailing stratiform precipitation is 
primarily to the south of the leading edge. Northwest 
of the radar site, convective precipitation between 
Bathurst and Melville Islands is also observed. The 
white circle in the precipitation images indicates a 
data-void zone in the radar volume scan. 

3 . SSM/I SIMULATIONS 

3.1 Radiative Transfer Calculations 

The liquid and ice precipitation water 
contents derived from the radar are converted to values 
of extinction coefficient, single -scatter albedo, and 
asymmetry factor at the SSM/I frequencies (19.35, 
22.235, 37, and 85.5 GHz) using the 

parameterizations in Kummerow and Weinman (1988). 
Absorption coefficients for molecular oxygen and 
water vapor are computed using the formulae of Liebe 
(1985). Microwave absorption by non-precipitating 
clouds is neglected in the present development. A 
surface skin temperature of 300 K and emissivities 
characteristic of land and ocean surfaces, where 
appropriate, are assumed. 

The upwelling microwave radiances at the 
SSM/I frequencies in two orthogonal planes of 
polarization at an incidence angle of 53.1° are 
computed using Eddington's second approximation 
along individual ray paths through the radar volume. 
The resulting field of microwave brightness 
temperatures at 2 km resolution at 37 GHz and 85.5 
GHz in the horizontal polarization are shown in Figs. 
2a and 2b, respectively. 

3.2 Sensor Response 

The spatial resolution of the SSM/I varies 
from about 55 km at 19.35 GHz to 14 km at 85.5 GHz, 
and antenna temperatures are sampled at intervals of 
12.5 km at 85 GHz, and at 25 km at all other 
frequencies. In order to simulate SSM/I resolution 
effects, the 2 km resolution images in Figs. 2a and 2b 
are convolved by the SSM/I 37 GHz and 85.5 GHz 
antenna patterns. The resulting images are then 
resampled at the SSM/I measurement locations. The 
final simulated antenna temperature images at 37 GHz 
and 85.5 GHz afe presented in Figs. 2c and 2d, 
respectively. Coincident antenna temperature 
measurements from the SSM/I are shown in Figs. 2e 
and 2f. 

4 . RESULTS AND DISCUSSION 

A general correspondence between 
precipitation regions in the simulated and observed 
SSM/I imagery is evident, and the magnitudes of the 
antenna temperatures are fairly consistent. However, 


the simulated antenna temperatures in the convective 
regions tend to be higher than those observed, while 
minimum simulated antenna temperatures over the 
stratiform precipitation region tend to be lower than 
those observed. These differences are most obvious in 
the 85.5 GHz imagery, for which scattering by 
precipitation-sized ice particles, and the associated 
lowering of upwelling radiances, is more pronounced 
than at 37 GHz. Errors in the partitioning of the liquid 
and ice precipitation components in the simulations 
might account for these discrepancies. Also it may be 
noted that simulated antenna temperatures at 37 GHz 
over the ocean tend to be somewhat lower than those 
observed. The lack of nonprecipitating clouds in the 
simulations might account for less emission and lower 
antenna temperatures over the low-emissivity ocean 
background. The sensitivity of simulated antenna 
temperatures to modeling assumptions will be the 
subject of future study. 
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